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OSCILLATORY  SUBSONIC  POTENTIAL  FLOWS  AROUND  THREE-DIMENSIONAL 
BODIES  AND  ITS  APPLICATION  TO  THE  CALCULATION  OF  DYNAMIC 
STABILITY  DERIVATIVES  OF  THE  AIRCRAFT 


Liu  Qiangang,  Wu  Changlin  and  Jian  Then 
(Northwestern  Polytechnical  University) 

Abstract 

This  article  introduces  a  unified  method  for  processing 
the  oscillatory  subsonic  potential  flows  around  three-dimen¬ 
sional  bodies  of  various  configurations.  The  major  feature  of 
this  method  is  the  employment  of  the  finite  element  method 
to  directly  explain  the  use  of  the  integro-differential 
equation  for  the  velocity  potential  on  the  surface  of  the  body 
derived  from  the  Green  theorem  to  obtain  the  velocity  poten¬ 
tial  distribution  on  the  surface  of  the  body.  Later,  we  further 
used  the  finite  difference  method  for  the  differential  of  the 
velocity  potential  so  as  to  obtain  the  pressure  distribution 
on  the  surface  of  the  body. 

Due  to  the  fact  that  theoretically  this  method  is  rela¬ 
tively  stringent,  when  used  for  the  calculation  of  flows 
around  bodies  with  complex  configurations  the  obtained  results 
were  relatively  accurate.  Because  of  this,  in  the  last  several 
years  its  application  has  become  more  and  more  widespread 
abroad.  Similar  basic  equations  from  related  reference  [3] 
were  utilized,  yet  there  were  dissimilarities  in  the  calcula¬ 
tion  of  the  aerodynamic  influence  coefficient.  This  article 
also  applies  this  method  for  the  calculation  of  dynamic  sta¬ 
bility  derivatives  of  the  aircraft  and  the  obtained  results 
are  in  agreement  with  the  experimental  results. 
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I.  Perturbation  Velocity  Potential  and  Boundary  Conditions 
If  we  take  a  right  angle  coordinate  system  such  as  the 
one  shown  in  fig.  1  (in  the  figure  the  OX  axis  direction  is 
identical  to  the  undisturbed  air  flow  speed  of  the  y  direc¬ 
tion)  under  small  perturbation,  there  is  the  following  per¬ 
turbation  velocity  potential  equation: 


Pig.  1, 


aHp  a*q >  a*«p  _L-/ rri^- 4-217 + 

7F-+  a?  +1F~- ^TV  U  +2U>*»t  +  <**  J 


In  the  formula, ¥  is  the  perturbation  velocity  potential,  Cqo 

is  the  speed  of  sound  in  the  undisturbed  air  flow  and  t  is 
the  time. 

If  we  can  use  the  following  equation  to  show  the  surface 
of  the  body: 

s  (x,  y,  z ,  t)  =0  (2) 

then  the  boundary  conditions  of  the  surface  of  the  body  is: 


v<p=>  - 


*(5- *«'-£) 


^  iw  '"ax; 

In  the  formula,  n  is  the  outer  normal  line  of  the  s. 

The  following  transformation  is  carried  out  for  the  var¬ 
iables  x,y,z,t  and  velocity  potential 

X-xftL,  Y-y/L,  Z-z/L 
T-Jit/L,  +-V/UL 


In  the  formula,  B=T/i-A£j//)t0  is  the  Mach  number  of  the  undis¬ 
turbed  air  flow  and  L  is  the  characteristic  length  of  the  body 


By  applying  the  Green  function  method  for  the  subsonic 
perturbation  flow  we  obtained  the  following  formula: 

«*♦(*•  i 


hXf{  CA*r^i)_Kr)  i 

Sw 


(5) 


In  the  formula,  S0  and  S separately  indicate  the  surface  of 
the  body  and  the  vortex  surface  of  the  body's  trailing  edge; 

21 4>  is  the  difference  of  the  velocity  potential  on  the  top  and 
bottom  of  the  vortex  surface;  p(X,Y,Z)  are  points  in  the  flow 
field;  q(X^,y.,Z^)  are  points  on  the  surface  of  the  body  or 
vortex  surface;  X^,Y^,Z^  are  integral  variables;  Nf)  is  the 
outer  normal  line  of  point  q;  R  is  the  distance  between  two 
points  of  p  and  q. 


(f  .“kTh  iZx-zy  (e) 


r  is  the  required  time  that  the  perturbation  is  emitted  from 
point  q  and  is  transmitted  to  point  p. 

t  +  (7) 

When  the  functions  indicated  in  the  brackets  takes  T^T-t;  E 
is  the  parameter  and  when  point  p  is  above  S0,E=l/2.  when 
point  p  is  outside  S_  and  not  above  S„  then  E=1 . 

D  W 


When  point  p  is  above  S  ,  formula  (5)  takes  <p  as  the 

□ 

integro-differential  equation  of  the  unknown  function.  It  can 
be  used  to  calculate  the  velocity  potential  distribution  of 
the  surface  of  the  body. 


It  is  worthy  to  note  that  in  formula  (5)  the  velocity 
potential  difference  on  the  top  and  bottom  of  the  vortex 
surface  is  actually  not  another  unknown  quantity.  It  can  be 
determined  by  the  difference  on  the  top  and  bottom  of  the 
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trailing  edge  of  the  surface  of  the  body: 

AM  « ,  7\)  -  A*(9«,  Ti—V’i)  (  8  > 

In  the  formula,  qTE  and  q  separately  indicate  points  on  the 
trailing  edge  of  the  surface  of  the  body  and  the  vortex  sur¬ 
face.  Moreover,  when  they  are  on  the  same  vortex  line?  M t 
is  the  required  time  of  the  vortex  point  flowing  from  point 

qTE  to  point  q: 

(ft) 

In  the  formula,  X  is  the  X  coordinate  of  point  q__ 

111  IE. 


The  above  integro-dif f erential  equation  can  use  the  finite 
element  method  to  seek  a  solution.  In  order  to  solve  this 
equation,  we  divided  the  surface  of  the  body  into  Ng  quadri¬ 
lateral  elements  and  took  the  &X  ftr,  ffl  in  each  of 
the  elements  to  be  constants  and  be  equal  to  their  values 
in  the  figure  centers  of  each  of  the  elements.  At  the  same 
time,  all  of  the  vortex  lines  are  divided  on  the  vortex  surface 
and  each  of  the  vortex  lines  and  series  of  elements  on  the  sur¬ 
face  of  the  body  are  joined.  The  length  of  the  vortex  line  can 

i6l 

take  any  effective  length.  Each  vortex  line  is  further  div¬ 


ided  into  many  quadrilateral  elements  and  the  value  of 
(Atf.  30  in  each  element  are  also  assumed  to  be  constants 
and  be  equal  to  their  values  in  the  figure  centers  of  each  of 
the  elements.  At  the  same  time,  we  also  considered  the  Kutta 
conditions  and  the  trailing  edge's  44*  of  the  surface  of  the 
body  could  be  considered  to  be  approximately  equal  to  the  4^ 
value  of  the  trailing  edge's  joined  element  figure  centers. 
Therefore,  formula  (8)  can  be  written  as: 


A4>*(r-T,»)=  (10) 

/ 

In  .the  formula,  i  and  j=l  and  2***Ng  is  the  serial  number  of  the 

element  on  the  surface  of  the  body;  h=l  and  2-**NTT  is  the 

w 
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serial  number  of  the  element  on  the  vortex  surface  (N  is  the 

w 

total  of  elements  on  the  vortex  surface);  K ..  is  the  coeffi- 

jn 

cient,  and  when  the  element  j  is  joined  with  the  trailing 
edge  of  the  surface  of  the  body  and  the  element  h  is  on  the 
trailing  edge  vortex  line  of  the  element  j,  Kjh=^l  (when  the 
element  j  is  on  the  top  surface  it  takes  a  and  when  on 
the  bottom  surface  it  takes  a  otherwise  Kj^=0?  A4> 

indicates  the  AP  of  the  h  element's  figure  center  on  the 
vortex  surface;  ^>j  indicates  the  p  of  the  j  element's  figure 
center  on  the  surface  of  the  body;tT;^is  the  time  required  for 
the  perturbation  to  be  transmitted  from  the  h  element ' s 
figure  center  on  the  vortex  surface  to  the  i  element’s  figure 
center  on  the  surface  of  the  body;  is  the  time  required  for 
the  vortex  point  to  flow  from  thesis  point  to  the  h  elements 
figure  center-  In  considering  the  use  of  the  above  approxima¬ 
tion,  when  calculating  4^,  the  trailing  edge's  coordinate  XTE 

of  the  surface  of  the  body  and  the  trailing  edge's  joined  X 
coordinates  of  the  surface  body  element  figure  center  can 
substitute  so  that  the  in  formula  (10)  can  be  written  as: 

(11) 

In  the  formula,  X^  and  separately  indicate  the  h  element 

on  the  vortex  surface  and  the  X  coordinate  of  the  j  element 
figure  center  on  the  surface  of  the  body. 

After  undergoing  the  above  process,  when  the  p  point  was 

on  the  SB,  formula  (5)  changed  into  a  constant  differential 
equation  with  a  constant  coefficient  and  linear  time  dif¬ 
ference  variable: 

+i<  r )  -  2  < c t  - + 2  c"+»c  t  -  *«)  +  s  D'*>( T  -  T"> 

i  i  i 

+  ^  2 Kn,Cik(.T  —  ■*u"”4*) 
i  h 

+  22  Ki»D,4,(  T  - Tu- 4«)  i  -  1 ,  2  - N . 

/  * 


(12) 


In  the  formula,  T;j  is  the  time  required  for  the  perturbation 
to  be  transmittd  from  the  j  element's  figure  center  of  the  sur¬ 
face  of  the  body  to  the  i  element's  figure  center?  fa  — 

A.\  J  '3T 

-JjijJ  )  L'j)  Cij,  tyjjQA'l  arS  the  aerodYnamic  influence 
coefficients: 


In  the  formula,  indicates  the  surface  of  the  j  element  on 

the  surface  of  the  body;  Sw^  indicates  the  surface  of  the  h 

element  on  the  vortex  surface;  b^  and  C  can  use  analytic 

formula  calculation  -  the  formula  for  calculating  is  the 

same  as  for  Cih  only  it  is  necessary  to  change  the  j  into  an 

h  in  the  formula.  For  approximate  calculation,  we  could  also 

use  D.  =R.  .C.  ,  and  D.,=R.,C..  .  In  these  formulas  R.  .  and  R., 
lj  ij  ij-  in  ih  xh.  ij  lh 

are  separately  the  j  element  on  the  surface  of  the  body  and 
distance  from  the'h  element's  figure  center  on  the  vortex 
surface  to  the  i  element's  figure  center  on  the  surface  of 
the  body. 

Equation  (12)  can  be  used  to  calculate  the  indeterminate 
constant  flow  around  bodies. 

II.  Calculation  Formula  for  Harmonic  Oscillatory  Flows  Around 
Bodies 
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When  a  body  revolves  around  its  fixed  constant  reference 
state  and  creates  harmonic  oscillation,  its  additional  per¬ 
turbation  velocity  potential  can  be  shown  as: 

4>(,Y,  Y,  Z,  T)=$(X,  Y,  ZVfl<Tt"-Jt,  (14) 

In  the  formula: 

Q=<bL/o-0=A'M-/P  (15) 

In  this  formula,  a  is  the  oscillatio  n  frequency;  K  is  the  con¬ 
version  frequency;  andjO  is  the  conversion  frequency  computed 
in  the  compressibility  effect. 

After  the  surface  of  the  body  is  divided  into  a  finite 
element,  formula  (14)  can  be  written  as: 

(16) 

When  formula  (16  is  substituted  into  formula  (12) ,  after  arrange¬ 
ment,  we  obtained: 

C&i;-C U,)=CbiO  <SJW,>  (17) 

In  the  formula,  §  „  is  the  Kronecker  function,  when  iifj 

then S  • .=0  ,  and  when  i=i  then  £ij=l. 

ij 


£//-«*'“"<  1  +:QR„)C„ 

(18) 

bt;=e'la*<ibt, 

(19) 

=  J]  *«#*'“**•***»<-“»>'*(  1  +*ORniCu 

(20) 

k 


When  the  body  creates  harmonic  oscillation,  the  equation 
for  the  surface  of  the  body  can  be  written  as  : 

*(.X,  Y,  Z,  T)=  Z-Z.AX,  Y)  *z(X,  Y)^  (21) 

In  the  formula,  corresponds  to  the  upper  surface;  “If"  cor¬ 

responds  to  the  lower  surface;  footnote U indicates  the  upper 
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surface, l  indicates  the  lower  surface,  and  Z.  (x,  y)is  the 
model  of  vibration.  Therefore,  the  boundary  condition  of  the 
oscillitory  flow  around  the  surface  of  the  body  can  be 
written  as: 

vs)=Nz(tKz+-y  jx)e'ioy’x  (22) 

In  the  formula,  N  is  the  cosine  of  the  constant  reference 

z 

state  of  the  outer  normal  line  N  on  the  surface  of  the  body 
and  the  oZ  axis  included  angle. 

In  general,  for  indeterminate  motion,  the  pressure  coef¬ 
ficient  can  be  calculated  according  to  the  following  formula: 

„  /  8  a<|>  ,1  *L\  (22) 

C'“  "  2  \M  aT  P  *X  ) 

and  for  the  harmonic  oscillatory  flow  around  the  body: 

C,=  CfeiaT  (24) 

From  formulas  (23),  (24)  and  (14)  we  obtained: 

C,=>~1  e1***  -fxae‘KX/9)  (2">; 

When  calculating  the  numerical  value,  the  finite  differ¬ 
ence  method  could  be  used  for  the  differential  in  the  above 
formula. 

After  ascertaining  the  geometric  parameters  of  the  surface 
of  the  body  and  the  boundary  conditions,  we  could  explain 
equation  group  (17)  and  obtain  the  value  ( j=l , 2  * • *N^) of 

the  N  dispersion  points  on  the  surface  of  the  body.  Further, 

B 

we  could  also  obtain  the  pressure  distribution  from  formula 
(25)  . 

III.  Its  Application  to  the  Calculation  of  Dynamic  Stability 
Derivatives  of  the  Aircraft 

Equation  (17)  has  been  applied  for  the  calculation  of  a 
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great  many  types  of  oscillatory  flows  around  bodies.  It  can 
be  used  to  calculate  aircraft  fluttering,  the  effects  of 
gusts  and  the  indeterminate  aerodynamic  forces  when  the  air¬ 
craft  makes  other  various  oscillatory  movements.  In  this 
section,  we  will  present  calculations  for  the  pressure  dis¬ 
tribution  and  dynamic  stability  derivatives  when  an  aircraft 
creates  pitch  oscillation  as  concrete  applications  of  this 
method . 

If  after  the  aircraft  has  perturbation  at  a  certain  equil¬ 
ibrium  and  tbs  revolves  around  its  center  of  gravity  creating 
pitch  oscillation,  its  model  of  vibration  can  be  indicated  as: 

z, )p  (26) 

In  the  formula,  X  is  the  x  coordinate  of  the  aircraft's  center 

b 

of  gravity. 

Based  on  the  method  introduced  in  the  above  section,  we 
can  obtain,  at  this  time,  the  variable  pressure  coefficient 
Cf>o  of  the  aircraft '  s  surface  and  carrying  out  the  integral 
of  this  variable  pressure  coefficient  obtain: 

c^cl+icl  (27) 

Cm—Cm  +  iCm  (28) 

In  the  formula,  letters  R  and  I  separately  indicate  the  real 
part  and  imaginary  part.  At  this  time,  the  aircraft's  lift 
and  pitching-moment  coefficient  are  separately.' 

C^RLCcue1**)  (29) 

C.-i?£(c-«<0T)  (30> 

In  the  formula,  "RL"  indicates  the  real  part. 

On  the  other  hand,  after  the  aircraft  has  perturbation 
and  creates  pitch  oscillation,  its  lift  and  pitching-moment 
coefficient  can  be  indicated  as: 
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Ct-Ct.o+CLT^— J+Cl^  „  J 

c.-c„+cm0  )  +C.i( *f) 

•  _  a 

In  the  formula,  cl  —  v  ^  tt 

In  the  above  two  formulas,  the  small  amounts  above  the 
second  order  are  omitted. ^ 

When  the  aircraft  has  pitch  harmonic  oscillation,  the 
increment  of  its  attack  angle  and  pitch  angle  is: 

a  -  0  <33> 


From  formulas  (29) -(33)  we  obtain: 
Ct=a«CCL.+j/C(Ct;+Cl,i)3  (34) 

C--«.CC«+iA:(C^  +C»i)3  (35) 


In  comparing  the  imaginary  parts  of  formulas  (27)  and  (34) 
and  formulas  (28)  and  (35)  we  obtain: 


(36) 

(37) 


Therefore,  CT*  and  C  •  can  be  obtained  under  fixed  constant 
L©  me  qq 

hypothetical  condition^.  The  value  of  the  dissimilar  con¬ 
version  frequency  times  of  C_~  and  C  7  can  be  obtained  from 

La  ma 

formulas  (36)  and  (37)  . 


IV.  Calculation  Examples 

We  applied  the  above  method  in  a  TQ-6  electron  computer 
at  the  Optical  Machine  Institute  of  the  Chinese  Academy  of 
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Sciences  in  Xian  for  the  calculation  of  aircraft  wings, 
ellipsoids,  the  left  distribution  of  the  pitch  harmonic  oscil¬ 
lation  created  by  a  certain  aircraft  (aircraft  wing,  fuselage, 
tail  assembly)  and  the  dynamic  stability  derivatives  CT-r  and 

C  x  of  an  aircraft.  The  results  were  as  follows  : 
ma 

1 .  Aircraft  Wings 

The  form  of  the  plane's  aircraft  wing  was  rectangular,  the 
aspect  ratio  was A  =2,  the  relative  thickness  was  c=0.001,  it 
revolved  on  a  x=c/2  axis  to  create  harmonic  oscillation  and 
when  in  M*,  =0  and  K=2  its  lift  distribution  was  as  shown  in 
figure  2. 


Fig.  2.  The  Lift  Distribution  of  the  Rectangular  Wing  Root 
Chord  Area  When  a  Rotating  x=c/2  Axis  Creates 
Harmonic  Oscillation 

1.  Results  of  this  article's  calculations 

2.  Reference  work  [2] 

It  can  be  seen  from  figure  2  that  the  results  of  this 
article's  calculations  are  identical  to  those  of  reference 
work  [2].  Because  the  divisions  into  pieces  were  relatively 
few,  the  lift  distribution  of  the  trailing  edge  section  near 
the  aircraft  wing  was  not  given  accurate  calculation  values  by 
reference  work  [2].  They  only  used  a  curve  to  indicate  its 
distribution  tendency.  Yet,  in  this  article,  because  we  used 
the  already  known  conditions  of  the  value  on  the  vortex 
surface,  although  the  division  into  pieces  were  few,  accurate 
calculation  results  were  still  obtained  for  the  trailing  edge 
section  near  the  aircraft  wing. 


2.  Ellipsoids 

The  ratio  of  the  major  axis  and  minor  axis  of  the  ellip¬ 
soid  was  8.  It  revolved  around  the  center  to  create  pitch 
harmonic  oscillation.  Its  pressure  distribution  is  shown  in 
chart  3 . 


Chart  3  The  Pressure  Distribution  When  an  a/b=8  Ellipsoid 
Creates  Pitch  Harmonic  Oscillation 

1.  Results  of  this  article's  calculations 

2.  Reference  work  [4  ] 

In  the  chart,  the  results  shown  by  the  symbol  "A  "  are 
the  analytically  interpreted  calculations  obtained  by  refer¬ 
ence  work  [4]  based  on  Helmholtz's  wave  motion  equation.  It 
can  be  seen  from  the  chart  that  given  the  situation  of  few 
divisions  into  pieces  (only  54  pieces  were  divided  on  half 
an  ellipsoid),  the  results  of  this  article's  calculations 
were  still  relatively  the  same  as  the  analytical  interpreta¬ 
tions. 


3.  Aircraft  Wings,  Fuselage  and  Tail  Assembly 
Below  we  will  give  the  test  results  of  a  certain  aircraft 
(aircraft  wings,  fuselage,  tail  assembly)  when  it  uses  dif¬ 
ferent  frequencies  to  revolve  around  its  center  of  gravity 
and  create  pitch  harmonic  oscillation.  The  outer  geometrical 
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form  of  the  aircraft  is  generally  as  is  shown  in  figure  4. 


F ig.  4 . 

Strictly  speaking,  the  dynamic  stability  derivative  should 
have  the  oscillation  frequency  be  close  to  zero.  Because  of 
this,  we  took  a  relatively  small  numerical  value  in  the  cal¬ 
culation  yet  this  did  not  cause  any  apparent  calculation 
errors  in  frequency  when  calculating  the  frequency.  Figures 
5,6  and  7  separately  show  the  lift  on  certain  sections  of 
the  aircraft  wings  and  tail  surface,  and  the  pressure  distri 
bution  on  a  certain  meridian  along  the  fuselage  when  the  air 
craft  used  these  frequencies  to  create  pitch  harmonic  oscil¬ 
lation. 


<«>  *> 


Fig.  5.  The  Lift  Distribution  on  a  Certain  Section  of  the 
Aircraft  Wing  When  a  Certain  Aircraft  Revolves 
Around  Its  Center  of  Gravity  and  Creates  Pitch 
Harmonic  Oscillation 
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Fig.  6.  The  Lift  Distribution  on  a  Certain  Section  of  a  Flat 
Wing  When  a  Certain  Aircraft  Revolves  Around  Its 
Center  of  Gravity  and  Creates  Pitch  Harmonic 
Oscillation 


ig.  7.  -The  Pressure  Distribution  on  a  Certain  Meridian  Along 
the  Fuselage  When  a  Certain  Aircraft  Revolves  Around 
Its  Center  of  Gravity  and  Creates  Pitch  Harmonic 
Oscillation 


See  table  1  for  the  calculation  results  of  the  aircraft's 
ci/«,  and  dynamic  stability  derivative  Cm&  when  M—  =0.6. 


K 

!  Cm  /  a8 

1  c  - 

»• 

0.005 

-  0.0545 

|  -  1.797 

0.01 

|  -0.109 

-  1.757 

Table  1 

—  (  1  ' 

In  the  calculations,  the  aircraft's  C  .  uses  -9,119v  ‘  . 

me 

It  can  be  seen  from  table  1  that  when  Mach  =  0.6,  the  air- 

00  9 

craft's  C  7  can  use  -  1.787. 
ma 

Therefore,  in  recently  published  works,  we  still  have  not 
seen  theoretical  calculations  or  experiment  results  for  the 
indeterminate  aerodynamic  force  of  aircraft  wings,  fuselage 
and  tail  assembly.  The  above  calculations  in  this  article  are 
still  not  suitable  material  to  offer  for  comparison.  Yet,  at 
present,  we  already  have  test  results  for  the  dynamic  stab¬ 
ility  derivatives  of  the  aircraft.  See  table  2  for  the  dynamic 
stability  derivative  test  values  of  the  aircraft  when  at 

Mach_  *0.6. 

00 


C  " 

!  c-* 

t 

-1.864 

-9.529 

(2)*»« 

-  1.787 

-9.119 

Table  2 

1.  Test  value 

2.  Calculated  value 


It  can  be  seen  from  table  2  that  the  test  values  of  the 
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aircraft's  dynamic  stability  derivatives  are ’quite  close  to 
the  calculated  values  of  reference  work  (1) .  This  also  ex¬ 
plains  in  one  sense,  this  article's  calculated  results  of 
the  indeterminate  aerodynamic  force  of  the  aircraft  wings, 
fuselage  and  tail  assembly. 

V  Conclusion 

This  article  introduced  a  unified  method  for  processing 
the  oscillatory  subsonic  potential  flows  around  three  dimen¬ 
sional  bodies  of  various  configuration.  It  also  presented 
calculation  results  for  the  indeterminate  aerodynamic  force 
and  dynamic  stability  derivatives  of  a  single  wing,  ellip¬ 
soids,  aircraft  wings  and  tail  assemblies  when  different 
frequencies  were  used  to  create  pitch  harmonic  oscillation. 
Actual  calculations  showed  that  the  application  range  of 
this  method  is  quite  wide  and  can  be  used  for  the  calculation 
of  the  oscillatory  flows  around  bodies  of  various  different 
and  complex  configurations.  Moreover,  the  amount  of  calcula¬ 
tion  work  is  relatively  small  and  the  results  are  relatively 
accurate.  At  present,  then,  this  is  a  good  method  for  cal¬ 
culating  oscillatory  subsonic  flows. 
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OSCILLATORY  SUBSONIC  POTENTIAL  FLOWS 
AROUND  THREE-DIMENSIONAL  BODIES  AND  ITS 
APPLICATION  TO  THE  CALCULATION  OF  DYNAMIC 
STABILITY  DERIVATIVES  OF  THE  AIRCRAFT 

Liu  Qiangang,  Wu  Changlin  and  Jian  Zheng 
( Northwestern  Polyiechmcal  University) 

Abstract 

A  general  formulation  for  oscillatory  subsonic  potential  flovrs  around  three- 
dimensional  bodies  of  various  configuration  and  its  application  to  the  cal¬ 
culation  of  dynamic  stability  derivatives  of  the  aircraft  are  presented.  By 
applying  the  Green  function  method,  we  obtained  an  integro-differential  equa¬ 
tion  relating  the  perturbation  velocity  potential  to  its  normal  derivatives  on 
the  surface  of  the  body.  In  order  to  solve  this  equation,  the  surface  of  the 
body  and  its  wave  are  divided  into  small  quadrilateral  elements.  The  unknown 
<t>  and  its  derivatives  are  assumed  to  be  constant  within  each  element.  Thus  the 
integro-differential  equation  reduces  to  a  set  of  differential-delay  equations  in 
time.  This  set  of  equations  can  be  used  as  the  basis  of  a  general  method  for 
the  fully  unsteady  flow  calculation.  For  oscillatory  subsonic  potential  flow,  this 
set  of  equations  further  reduces  to  a  set  of  linear  algebraic  equations  which  is 
solved  numerically  to  yield  the  values  of  <j>/  at  the  centroid  of  each  element. 
The  pressure  coefficient  is  evaluated  by  the  finite  difference  method.  The  lift 
and  the  moment  coefficients  are  determined  by  numerical  integration  of  the 
pressure  coefficient.  The  dynamic  stability  derivatives  are  obtained  from  the 
imaginary  parts  of  the  lift  and  the  moment  coefficients. 

The  formulations  in  this  paper  are  embedded  into  a  general  computer  pro¬ 
gram.  Several  typical  numerical  results  have  been  obtained  by  means  of  this 
program.  Figure  2  shows  the  distribution  of  lift  coefficient  CL  along  the  middle 
section  for  a  rectangular  wing  oscillating  in  pitch  with  X  “2,  t  *0.001,  M-=*0, 
K’"  2  .The  result  is  identical  to  the  original  calculation  by  Morinoc,J.  Figure  3 
shows  the  distribution  of  pressure  coefficient  Cr  for  a  harmonically  oscillating 
spheroid  with-|-=8,  M«  =  0.5,  K  —  2.  The  result  is  in  good  agreement 
with  the  analytical  solution  of  wave  equation143. Figures  5 ,  6,  7  show  the  di¬ 
stributions  of  lift  coefficient  at  various  stations  of  an  aircraft  (wing-body- 
tail  combination) oscillating  in  pitch  with  M.  =  0.6,  K  =  0.005, 0.01.  Table  2  shows 
the  dynamic  stability  derivatives  C;.«,  C„;  of  the  aircraft.  The  results  are 
in  good  agreement  with  the  experimental  data. 
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AN  INTERPOLATION  MIXED  DIFFERENTIAL  METHOD  FOR  TRANSONIC  LARGE 
DISTURBED  SYMMETRICAL  POTENTIAL  FLOW  AROUND  AIRFOIL 

by  Ling  Heyao 

(Design  Department  of  Hongan  Aircraft  Company) 

Abstract 

By  extending  an  interpolation  mixed  differential  methodri] 
for  transonic  small  disturbed  steady  potential  flow  to  the  tran¬ 
sonic  large  disturbed  steady  potential  flow,  we  proposed  an 
interpolation  mixed  differential  method  for  solution  of  the  ex¬ 
act  equation  for  transonic  potential  flows  in  the  local  speed 
coordinate  system-  In  numerical  illustration  of  this  method,  the 
pressure  distributions  of  the  double  arc  airfoil  and  NACA0015 
airfoil  in  symmetrical  state  are  computed  and  compared  with  the 
data  of  the  experiments  [2,3]  and  the  results  of  the  computa¬ 
tion  for  the  double  arc  airfoil  by  the  small  disturbed  mixed 
differential  method.  The  results  were  close.  The  computations 
proved  that  the  interpolation  mixed  differential  scheme  is 
stable  and  convergent.  This  paper  has  solved  the  difficulty  of 
computing  the  Mach  w  number  which  near  1. 

I.  Preface 

Since  1970  when  Murman  and  Cole[4]  proposed  the  use  of  an 
interpolation  mixed  differential  method  for  solving  the  tran¬ 
sonic  small  disturbed  potential  flow  equation,  various  different 
methods  have  continued  to  appear  for  the  solution  of  transonic 
flow  and  development  has  been  very  fast.  Beginning  from  the  end 
of  1973,  using  reference  [4]  as  a  basis.  Professor  Luo  Shijun 
further  developed  and  extended  the  computation  and  application 
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of  transonic  steady  disturbed  potential  flow  and  obtained  a 
great  many  resultsfl].  In  order  to  extend  the  application  range 
of  the  mixed  differential  method,  reference  [6]  extended  the 
mixed  differential  scheme  for  the  transonic  small  disturbed 
potential  flow  to  the  transonic  large  disturbed  potential  flow. 
The  revolving  mixed  differential  scheme  was  proposed  whereupon 
computation  results  were  obtained  for  the  obtuse  flow  around 
the  body  yet  the  amount  of  computation  work  was  relatively 
large.  Because  of  this,  this  article  proposes  a  new  computation 
method  which  is  an  interpolation  mixed  differential  method  for 
the  local  speed  coordinate  system. 

II.  Speed  Potential  Equation  in  the  Local  Speed  Coordinate 
System 

In  a  plane  flow,  the  exact  speed  potential  equation  is: 

(a* — u1)  <|>„  +  (a1 — ul)  <j>„  —  2uvfy,y  =  0  (  j  ) 

In  the  formula,  a  indicates  the  speed  of  sound;  u  and  v 
separately  indicate  the  speed  component  along  the  x  and  y 
coordinate  axes;  and  <$)  (x,  y)  indicates  the  speed  potential. 

In  the  son  local  speed  coordinate  system,  s  is  the  tangent 
of  the  0  point  on  the  flow  line  and  n  is  the  normal  line  of  the 
0  point  on  the  flow  line  (see  fig.  1) .  In  the  local  speed 
coordinate  system,  the  local  speed  q  is  identical  to  the  s  axis. 
Because  of  this,  there  is  no  division  of  speed  in  the  n  axis 
direction. 
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Fig.  1.  The  Local  Speed  Coordinate  System 
Key:  1.  Flow  line 

Based  on  this,  when  equation  (1)  is  converted  to  the  local  speed 
coordinate  system,  it  possesses  the  following  equation  form: 

(a1  —  gl)4>„  +  a1<t>,.=  0  (  2  ) 

2 

When  the  two  sides  of  formula  (2)  are  divided  by  a  ,  we  obtain: 

( l -M2)4>«+4>,...=  o  (3) 

In  the  formula,  M=q/a  is  the  local  Mach  number. 

Equation  (3)  is  a  non-linear  partial  differential  equation. 
2 

Given  that  the  1-M  item  can  be  positive, negative  or  zero,  it 
is  also  an  equation  which  can  be  an  equational  form  in  a  flow 
line  for  an  elliptical  form,  hyperbolic  form  or  parabolic 
form.  Because  of  this,  equation  (3)  is  a  mixed  form  and  the 
solution  of  the  flow  field  is  a  mixed  problem. 

In  the  supersonic  or  sonic  range,  equation  (3)  is  a  hyper¬ 
bolic  form  or  parabolic  form  with  an  existing  group  of  char¬ 
acteristic  lines.  In  the  local  coordinate  system,  the  charact¬ 
eristic  line  is  symmetrical  to  the  local  speed  axis  and  further¬ 
more  the  characteristic  line  makes  a  local  Mach  angle  H  with 
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the  included  angle  of  the  local  speed.  Thus,  the  character¬ 
istic  line  is  exactly  the  local  Mach  angle.  In  the  small  dis¬ 
turbed  plane  potential  flow,  the  characteristic  line  is  sym¬ 
metrical  to  the  basic  coordinate  axis,  and  the  two  situations  of 
the  characteristic  line  equation  forms  in  each  coordinate 
system  are  completely  identical.  Because  of  this,  in  the  same 
way  we  conclude  that  the  stability  of  the  differential  scheme 
in  the  small  disturbed  potential  flow  is  appropriate  for  the 
large  disturbed  plane  potential  flow  in  the  local  speed  coordin¬ 
ate  system. 

III.  The  Fitting  and  Interpolation  of  the  Speed  Potential 
Surface 

The  speed  potential  surface  in  the  binary  flow  is  a  complex 
surface.  It  is  necessary  to  use  a  single  function  to  express 
the  formula  exactly  and  this  is  actually  impossible.  We  should 
use,  based  on  the  network's  nodal  point  distribution,  any 
"speed  potential  surface  piece"  to  describe  it.  Within  each 
speed  potential  surface  piece  the  smoothness  requirements  in 
the  speed  potential  curved  surface  piece  should  be  guaranteed. 

In  this  way,  after  joining  with  a  suitable  form  of  a  speed 
potential  surface  piece,  nearly  any  form  of  speed  potential 
surface  piece  can  express  it  and  thus  the  required  accuracy  can 
be  attained. 

Based  on  reference  work  [7],  we  can  find  the  boundary  line 
for  the  ruled  surface  interpolation  formula  of  the  line.  The 
angular  point  of  the  already  known  ruled  surface  is  00,01,10,11 
(see  figure  2).  We  can  form  a  set  of  boundary  lines: 

o«  =coo  our  f„  ' 

If,  . 

lK-oo  wr  f„  ' 

If,  . 


(4) 

(5) 
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Fig.  2.  The  Position  Information  of  the  Ruled  Surface  Angular 
Points 

After  attaining  these  two  lines,  the  identical  parameter 
value  n  points  on  their  surface  are  viewed  as  corresponding 
points  and  lines  are  drawn  in  each  of  the  pairs  of  correspond¬ 
ing  points  in  the  same  manner.  When  this  type  of  line  becomes 
the  m  line  on  the  ruled  surface,  we  obtain  the  m  line 
formula  : 


mn-CF,  FJ 


0  n 

U» 


(6) 


When  we  make  this  type  of  m  line  the  "generating  line"  and 
the  slip  on  the  two  "base  lines"  of  formulas  (4)  and  (5) 
causes  the  n  to  change  from  0  to  1 ,  we  obtain  the  total  m  line 
of  the  ruled  surface.  This  is  also  the  scanning  of  the  entire 
ruled  surface.  Thus,  formula  (6)  is  an  equation  the  ruled 
surface. 

Using  formulas  (4)  and  (5),  we  can  write  formula  (6)  as: 


(rrm)  —  CF, 


F,3f  00 
.  10 


(7) 


The  entire  element  in  the  divalent  square  matrix  on  the 
right  side  of  the  formula  is  the  constant  vector  and  they  are 


all  position  information  of  the  angular  points  on  the  ruled 
surface  angular  points.  Therefore,  this  square  matrix  is 
called  the  angular  point  information  square  matrix.  The  inform¬ 
ation  in  the  square  matrix  is  easily  provided  by  actual  pro¬ 
blems  and  formula  (7)  is  also  convenient  for  numerical  com¬ 
putations. 


Fig.  3.  Double  Cubic  Surface 

If  the  known  surface  boundary  lines  are  cubic  functions, 
then  they  accordingly  form  a  ruled  surface  equation  which  in 
turn  can  form  a  double  cubic  surface  interpolation  formula 
(see  fig.  3) .  Without  undergoing  deductions,  an  equation  is 
given  similar  to  formula  (7) : 


(mn)«CF,  Fj  Gt  (?,3 

00  01  00.  01. 

F, 

10  11  10.  11. 

F, 

O 

o 

1 

o 

*-* 

1 

o 

o 

1 

■ 

o 

»-* 

• 

• 

<?. 

\Qm  Hal  10m»  11<»* 

Functions  F  F,G  and  G,  in  formulas  (7)  and  (8)  are  called 
o  l  o  I 

mixed  functions f 7 ] .  The  complete  elements  in  the  tetravalent 
square  matrix  of  formula  (8)  are  constant  vectors.  They  are  all 
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position  vector,  tangent  vector  and  torsional  vector  informa¬ 
tion  on  the  angular  points. 

The  tangent  vectors  and  torsional  vectors  on  the  surface’s 
(mn)  four  angular  points  can  use  the  differential  method  com¬ 
putation  (see  fig.  3) .  The  tangent  vector  differential  form¬ 
ula  on  the  ij  nodal  point  is: 


f  »_(  «  +  1  ■  J )-( » -  1 »  ■?)  -(-  o(A) 

(i,  /  +  0-0  »  L  ~li-+  o  (a) 


ij 


(9) 

(10) 


Z: m  and  An  are  the  nodal  points  of  the  two  parametric  change 
directions . 

The  torsion  vector  differential  formula  on  the  ij  nodal 
point  is: 


r( »  + 1 ,  /  —  i)—(  >•  - 1 ,  /  - 1 ) 

l  A»l|.lfM 

iiiijt.d) 


+  An// L  A»«f. +  Am<p;., 

(  i  +  1 1  ;  -f  1 )  —  ( i  —  l , 

+  A»»»„-«.i 


(11) 


The  tangent  vectors  and  torsion  vectors  on  the  other  three 
angular  points  can  be  computed  in  a  similar  manner. 

The  differential  formula  for  the  tangent  vector  and  torsion 
vector  on  the  boundary  can  be  formed  according  to  the  imbed¬ 
ding  boundary  condition  and  exterior  extent  topological  methods 
mentioned  in  reference  [  1 ] .  We  will  not  go  into  greater  detail 
here . 
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The  interpolation  accuracy  of  interpolation  formula  (8) 
is  higher  than  that  of  formula  (7) .  Formula  (7)  only  guarantees 
the  continuance  of  the  positions  on  the  adjacent  surface 
boundaries;  the  boundary  slope  is  not  continuous  and  shows  a 
surface  bending;  the  adjacent  curved  surface  boundary  line  is 
also  not  continuous  in  the  tangent  vector  on  the  angular  point 
and  shows  a  curved  surface  inflection  thus  forming  a  space 
broken  line.  However,  formula  (8)  not  only  guarantees  the  con¬ 
tinuance  of  the  position  on  the  adjacent  curved  surface  bound¬ 
ary  but  also  ensures  the  tangent  vector  continuance  of  the 
boundary  slope  and  boundary  line  on  the  angular  points. 

IV.  The  Differential  Scheme  and  Boundary  Conditions 

The  differential  scheme  explained  by  the  mixed  differen¬ 
tial  of  a  non-linear  mixed  partial  differential  equation  must 
be  based  on  the  flow  field  at  subsonic,  sonic  and  supersonic 
speeds.  Then  equation  (3)  can  separately  choose  the  different 
patterns  of  an  ellipsoid,  paraboloid  and  hyperbolic.  In  order 
for  the  differential  scheme  of  the  linear  hyperbolic  type 
equation  to  satisfy  the  Courant-Friedrichs-Lewy  stability  con¬ 
ditions  in  the  local  speed  coordinate  system,  the  dependent 
area  of  the  difference  equation  must  be  greater  than  the 
dependent  area  of  the  differential  equation.  If  we  use  the 
center  difference  on  the  n  axis  and  use  the  upstream  one  side 
difference  on  the  s  axis,  then  the  dependent  area  of  the 
difference  equation  is  necessarily  greater  than  or  equal  to 
the  dependent  area  of  the  differential  equation.  Also,  the 
stability  conditions  are  always  satisfied.  The  perturbation  in 
the  subsonic  flow  field  can  be  disseminated  in  every  direction 
and  conversely  it  is  also  dependent  on  the  points  all  around 
the  flow  field.  Therefore,  both  the  n  axis  and  s  axis  use  the 
center  difference  scheme  and  are  able  to  exactly  describe  the 
physical  characteristics  of  subsonic  perturbation. 


Taking  the  ij  nodal  point  as  the  basic  point,  they  sep¬ 
arately  interpolate  the  step  length  on  the  local  speed  coord¬ 
inate  system  as  the  corresponding  speed  potential  values  of 

£  i ,  fi2'  ^3  and  $4  of  ^  ' 

The  center  difference  form  is  used  for  the  tangent  flow 
field  (see  fig.  4). 

(l2) 

The  upstream  one  side  difference  form  is  used  for  the  inner 
points  of  the  supersonic  tangent  flow  field  (see  chart  4) . 


The  center  difference  form  is  always  used  for  the  normal 
flow  field  (see  fig.  4) . 


0  (Al) 


^  <2) 


Fig.  4. 
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Fig.  4.  The  Difference  Scheme  for  the  Inner  Points  of  the  Flow 
Field 

Key:  1.  Center  form 

2.  One  side  form 


l* 
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Fig.  5.  The  Difference  Scheme  on  the  Surface  of  the  Body 

Key:  1.  Center  formula 

2.  One  side  formula 

The  boundary  condition  for  the  surface  of  the  body  in 
formula  (3)  is: 

0  05) 

Because  the  boundary  line  for  the  surface  of  the  body  is  a 
flow  line,  the  speed  on  the  normal  surface  of  the  body  is 
naturally  0  and  the  local  speed  q  must  be  in  contact  tangent 
with  the  surface  of  the  body.  When  the  boundary  conditions 
of  the  surface  of  the  body  are  imbedded  into  the  normal  diva¬ 
lent  partial  derivative  we  obtain  the  difference  formula  on 
the  surface  of  the  body: 

|r (<!>«- 4>n)+ 0(4)  (l',} 
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The  tangent  divalent  partial  derivatives  separately  use 
formulas  (12)  or  (13)  based  on  whether  they  are  subsonic  or 
supersonic . 

When  this  article  concretely  processed  the  boundary  but 
had  not  yet  satisfactorily  obtained  the  actual  surface  of 
the  body,  the  position  was  shifted  on  to  the  chord  line.  The 
tangent  direction  was  still  identical  to  the  tangent  direction 
on  the  surface  of  the  body.  That  is  vector  s  on  the  chord  line 
was  parallel  to  vector  s1  on  the  surface  of  the  body  (see 
Fig.  6) . 


Fig.  6.  Simplified  Schematic  of  the  Boundary 

The  local  speed  q  in  the  Mach  number  square  represented 
form  always  uses  the  center  difference  scheme: 

q  „  A"**-*  0  (A*)  07) 

In  the  flow  field,  each  point  of  the  local  speed  q  and  the 
basic  coordinate  axis  included  angle  &  is  different.  Based  on 
speed  components  u  and  v  in  the  basic  coordinate  system  we  can 
obtain: 


6  *tg" 


v 

u 


(18) 


In  the  formula,  both  u  and  v  use  the  center  difference  form: 


u 


h+j wzALluL  +  0  (  A) 

A*i-i  + Axi 


v 


by  1.1  + by  i 


+  0(A) 


(19) 


V.  The  Difference  Equation  and  Its  Analytical  Method 

Given  that  the  flow  field  is  subsonic  or  sonic  and  super¬ 
sonic,  separating  composite  formulas  (12)  —  (17) ,  we  obtained 
the  mixed  difference  equation  of  the  flow  field. 


For  the  inner  points  and  points  on  the  surface  of  the  body 

2 

in  the  subsonic  flow  field,  1-M  ^O. 


(  1  —  M*)(4»i  —  2<|>i/+<frt)  + 


j2(<f»* — 4*/#)  "  o 

l(<i>3—  2$»f+$«)  =“ 


0 


(l) 

rtA  (2) 


(20) 


Key:  1.  Points  on  the  surface  of  the  body 

2.  Inner  points 

For  the  inner  points  and  points  on  the  surface  of  the  body 
in  the  supersonic  flow  field  (including  the  sonic  points) , 

1-M2  ^  O 


( 1  —  M‘)  (♦» — 2$i + 4 + 


2(<K -$*/)■■  0 

.(♦»  —  24>f/  +  $4) 


0 


(1) 

(21) 

(2) 


Key:  1.  Points  on  the  surface  of  the  body 

2.  Inner  points 
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The  equation  of  sonic  discriminate  difference  is: 


1  —  M*= 


r  4A*5i 

l  (4»s— 4».)^ 

L  2  +  Mi  /  2  J 

In  the  formula,  q*,  =1  and y =1.4. 


(22) 


Besides  i>  i j ,  difference  equations  (20)  and  (21)  were 
both  obtained  by  previous  field  speed  potential  interpolation. 
Because  of  this,  we  used  the  point  relaxation  simple  iteration 
solution.  Its  formula  is  differentiated  as: 


For  1-M  >  0, 


24>J 

(<l»3  +  <)>,))  f*I.&  (2) 


( 2;: ) 


Key:  1.  Points  on  the  surface  of  the  body 

2.  Inner  points 


For  1-M2  0, 


-  1 


+  M 


TC(1  -M*)(4. 


i  *~2<|>j)+| 


2<i>«) 

(<l>j  +  <i><))  (2) 


(24) 


Key:  1.  Points  on  the  surface  of  the  body 

2.  Inner  points 

In  order  to  quicken  the  computation  convergence  speed, 
the  relaxation  operation  formula  was  applied  in  the  relaxation 
iteration  computation  process: 

4. 17'  =oo,T,,(r +4>,‘7':'  (i  —  °> )  (25) 

In  the  formula,  is  the  iteration  computation  result  of  the  n 
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sequence  having  not  undergone  relaxation  processing,  $  ^  and 
$  v  '  are  the  iteration  computation  results  of  the  n  sequence 
and  n-1  sequence  after  having  gone  through  relaxation  process¬ 
ing;  W  is  the  relaxation  factor,  its  value  is  determined  by 
requirements  and  it  generally  uses  W  ^  but  this  paper's 
computations  generally  takes  it  as  1. 

The  computation  of  the  initial  field  uses  the  undisturbed 
homogeneous  flow  field  or  the  computation  result  field  near  the 

Mach  number.  The  formula  for  the  undisturbed  homogeneous  flow 
00 

field  is: 


The  distant  field  boundary  speed  potential  value  without  ex 
ception  uses  the  undisturbed  homogeneous  flow  field  value. 

The  pressure  coefficient  formula  is: 


-{[  1  -  (Y-*-  —  (q'-ql)  )  T  -  t  -  1  | 


In  the  formula,  local  speed  q  uses  the  center  difference  formula 

VI.  Numerical  Illustrations  and  Their  Analysis 

This  article  writes  on  two  numerical  illustrations.  One  is 
the  sharp  nosed  double  arc  airfoil  and  the  other  is  the  blunt 
nosed  NACA0015  airfoil.  Their  relative  thicknesses  are  0.06  and 
0.15  respectively. 


The  computation  was  done  on  the  TQ-6  computer  which  has 
close  to  a  1  million  time  operational  speed.  Two  types  of 


computation  networks  were  used:  the  42x25  and  the  82x40  were 
used  to  test  the  effects  of  different  network  densities  on  the 
computation  results  and  iteration  convergences.  The  computa¬ 
tions  showed  that  the  different  network  densities  had  no  effects 
on  computation  stability.  The  figure  results  of  these  computa¬ 
tions  all  used  the  42x25  network. 

The  computation  results  uniformly  used  two  similar  wing 
surfaces  with  pressure  coefficient  errors  of  Ap  ^  0.0001. 

See  figs.  7  and  8  for  the  computation  results. 


Fig.  7.  The  Pressure  Distribution  of  the  Double  Arc  Airfoil 

Key :  1 .  Our  method 

2.  Test  values  [  3  ] 

3.  Small  disturbed  difference  method 

Fig.  7  appends  the  test  values  and  small  disturbed  dif¬ 
ference  method  computation  values  from  reference  f 3 1 .  The 
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results  in  this  paper  coincide  very  well  with  the  test  values 
and  are  actually  closer  to  the  test  values  than  the  small  dis¬ 
turbed  difference  method. 


Fig.  8.  The  Pressure  Distribution  of  the  NACA0015  Airfoil 

Key:  1.  Our  method 

2.  Test  values  [2] 

Fig.  8  appends  the  test  values  from  reference  [ 2 1 .  When 
in  a  subcritical  state,  the  results  of  this  article  coincide 
very  well  with  the  test  values.  Yet,  in  a  supercritical  state, 
the  pressure  coefficient  computation  values  are  higher  than  the 
test  values.  When  at  a  supercritical  state,  the  airfoil’s 
boundary  layer  is  even  larger.  This  article  did  not  consider 
the  effects  of  the  increased  thickness  of  the  boundary  layer 
and  therefore  the  computations  tend  to  be  high. 
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Due  to  the  fact  that  in  the  computation  process  of  this 
paper  the  stability  was  very  stable  there  was  no  acute  oscil¬ 
lation  divergence.  Therefore  the  computation  initial  values 
were  uniform  for  the  flow  field. 

In  reference  [1]  the  computation  of  the  Mach  number  was 
close  to  1  and  therefore  the  oscillation  divergence  which  used 
the  computation  of  the  initial  field  was  the  result  field  of 
the  above  computed  Mach  M  number.  Moreover  the  computed  Mach^ 
number  increment  was  very  small  and  ^  Mm  =0.001.  That  is, 
with  the  use  of  these  types  of  strict  measures  the  computed 
Mach  M  number  only  reached  0.953.  Zheng  Youwen,  one  of  the  writers 
of  recently  composed  article  (1),  utilized  the  zero  initial 
field  as  well  as  the  control  relaxation  factor  method  and  com¬ 
puted  the  NACA0012  blunt  nosed  airfoil  whereupon  he  only 
reached  Mw  =0.95.  He  also  pointed  out  that  up  until  the  present, 
when  using  a  precise  potential  flow  equation  to  compute  the 
NACA001 2  airfoil,  after  MB  ^  0.9,  the  computations  were  diver¬ 
gent  and  no  convergent  computation  results  were  obtained. 
Nevertheless,  this  article  obtained  convergent  computation  re¬ 
sults  for  the  NACA0015  and  double  arc  airfoil  when  Mach  was 

00 

close  to  1  as  well  as  corresponding  Mach  M  numbers  of  0.999  and 
1.001  (the  results  for  the  double  arc  airfoil  are  not  shown  in 
chart  7)  .  Numerically,  the  two  Mactj,,  computation  results 
are  almost  identical  and  thus  we  have  overcome  the  difficulty 
of  computing  the  Mach  b  near  1. 

VII.  Conclusion 

This  method  is  applicable  for  subsonic,  transonic  and 
supersonic  flows. 

The  computation  accuracy  of  this  method  is  higher  than  the 
mixed  differential  method  for  transonic  small  disturbed  poten¬ 
tial  flow.  Because  the  equation  is  a  hypothesis  for  small  un¬ 
disturbed  perturbation,  with  the  expansion  of  the  computation’s 
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application  range,  we  can  compute  the  large  perturbation  flow 
field  of  a  blunt  nosed  body. 
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AN  INTERPOLATION  MIXED  DIFFERENTIAL 
METHOD  FOR  TRANSONIC  LARGE  DISTURBED 
SYMMETRICAL  POTENTIAI  FLOW  AROUND 

AIRFOIL 

Ling  Heyao 

(Design  Department  of  Hong  An  Aircraft  Company) 

Abstract 

By  extending  an  interpolation  mixed  differential  method  for^  transonic 
small  disturbed  steady  potential  flow  to  the  transonic  large  disturbed  steady 
potential  flow,  we  proposed  an  interpolation  mixed  differential  method  for 
solution  of  the  exact  equation  for  transonic  potential  flows  in  the  local 
speed  coodinate  system.  In  numerical  illustration  of  this  method  the  pressure 
distributions  of  double  arc  airfoil  and  NACA0015  airfoil  in  symmetrical 
state  are  computed  and  compared  with  the  data  of  experiments^’5-1  and 
the  results  of  the  computation  for  the  double  arc  airfoil  by  the  small 
disturbed  mixed  differential  method.  The  comparisons  are  shown  in  a  good 
approximation,  therefore  it  is  proved  that  the  interpolation  mixed  differential 
scheme  is  stable  and  convergent.  This  paper  has  found  a  way  out  of  the  diffi¬ 
culty  at  Mach  number  near  1, 
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THE  CALCULATION  OF  LIFT  AND  DRAG  CHARACTERISTICS  OF  SUBSONIC 
WINGS  WITH  WINGLETS 

by  Zhou  Renliang 

(Nanjing  Aeronautical  Institute) 

Abstract 

This  paper  uses  the  finite  fundamental  solution  to  divide 
the  spanwise  lattice  and  determine  the  spanwise  locations  of 
control  points  by  means  of  the  constant  roll-angle  method.  We 
calculated  the  lifts  of  rectangular  wings  with  different  winglets 
at  subsonic  speeds  and  calculated  the  induced  drags  by  using  the 
combined  flow  field  method.  From  the  calculations  of  various 
configurations  of  a  winglet,  we  have  found  out  some  rules 
affecting  the  lift  and  drag  characteristics  of  wings  and  picked 
out  a  favorable  configuration  from  them.  The  mechanism  of 
winglets  is  also  analyzed. 


Symbols 

A*  Degree  of  constant  roll-angle 

&  The  included  angle  of  the  roll-angle  ray  and  y  axis 

rf  Total  number  of  lines  of  the  lattice 

Total  number  of  rows  of  the  lattice 

5  Line  series  number  of  the  lattice 

6  Wing  chord  length 

/  Wing  semi span  length 
5  Winglet  angle 
a  Wing  attack  angle 
5  Wing  area 
A.  Wing  aspect  ratio 
dS  Lattice  area 
Y  Circulation  coefficient 
ACf>  Load  coefficient 
D  Normal  induced  speed  on  wing 
Vjp  Overtaking  flow  speed 
Cy  Wing  lift  coefficient 
C*i  Wing  induced  drag  coefficient 
^  Induced  drag  factor 
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I .  Preface 

A  winglet  is  an  aerodynamic  plan  with  a  simple  and  effect¬ 
ive  structure  as  well  as  a  new  technology  for  improved  aircraft 
performance.  The  winglet  which  is  fitted  on  the  wing  can  raise 
the  lift  coefficient  and  lower  the  induced  drag  coefficient  of 
a  subsonic  aircraft,  thus  raising  the  climbing  rate  and  econ¬ 
omizing  on  cruise  fuel  consumption.  Therefore,  at  present,  the 
winglet  has  already  been  given  serious  attention  in  aerodynamic 
research  and  by  aircraft  designers. 

Although  some  articles  and  test  results  have  been  published 
on  winglet  research  still  there  have  not  been  very  many.  Further¬ 
more,  there  has  not  been  ample  theoretical  research  and  aero¬ 
dynamic  mechanism  analysis  on  the  winglet. 

This  paper  which  employed  the  finite  fundamental  solution 
method  recommended  in  reference  [1],  the  constant  roll-angle 
method  to  divide  the  spanwise  lattice  and  to  determine  the 
spanwise  locations  of  control  points,  and  the  combined  flow  field 
method  to  calculate  the  induced  drag,  obtained  satisfactory  lift 
and  induced  drag  values.  On  this  foundation,  we  calculated  and 
compared  the  lift  and  drag  characteristics  for  wings  with  winglets 
of  various  configurations.  In  this  way,  we  found  out  some  rules 
affecting  the  aerodynamic  layout  of  the  winglet  and  picked  out 
a  favorable  configuration  from  them.  We  also  analyzed  the  aero¬ 
dynamic  mechanism  of  the  wing's  vortex  elimination  action  and 
reduction  of  reduced  drag  by  means  of  the  circulation  distribu¬ 
tion  chart  and  distribution  chart  for  the  circulation  change 
rate. 

II.  Computation  Methods 

1.  Calculation  of  the  wing  lift  coefficient 

We  used  the  constant  roll -angle  method  to  divide  the  span- 
wise  lattice  and  to  determine  the  spanwise  locations  of  control 
points.  The  principle  of  the  constant  roll-angle  method  is  as 
follows : 
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Taking  the  semispan  length  2*  as  the  radius,  the  center  of 
the  circle  is  placed  on  the  coordinate  origin  to  create  a  1/4 
circle  and  beginning  from  the  y  axis  each  AQ  angle  is  a  meri¬ 
dional  ray  and  thus  we  can  obtain  a  group  of  meridional  rays. 
The  &&  angle  is  determined  by  the  following  formula: 


AB  = 


(t  ) 


The  included  angle  of  each  ray  and  y  axis  is  ©: 


e=KAd  K  =  1,  2-2N  ( 2 ) 


When  the  point  of  intersection  of  these  rays  and  the  arc 
are  projected  on  the  y  axis,  K  is  the  projection  of  the  odd 
numbered  ray  point  of  intersection  on  the  y  axis.  This  becomes 
the  spanwise  location  of  control  points  and  K  is  the  even 
numbered  ray  point  of  intersection  on  the  y  axis  which  becomes 
the  lattice  dividing  line.  For  the  chordwise,  we  used  the 
partitioned  lattice  and  chordwise  location  of  control  points 
taken  on  the  1/2  chord  line  of  the  lattice.  On  the  half  wing,  the 
common  M  X  N  lattice  becomes  a  winglet  section  beginning  from 
a  certain  row  to  the  wing  tip  causing  it  to  have  the  required 
warp  angle  and  thus  forming  the  unflat  surface  of  a  wing  tip 
with  a  winglet. 

On  each  lattice  there  is  arranged  a  compressible  horse  shoe 
vortex,  there  is  added  a  chordwise  location  placed  on  the  ad¬ 
vancing  edge  of  the  lattice  and  two  free  vortex  separately  fol¬ 
low  the  overtaking  flow  from  the  two  ends  of  the  lattice's 
advancing  edge  and  extend  out  towards  infinity.  Each  control 
point  satisfies  the  flow  boundary  conditions.  It  can  solve  the 
linear  algebraic  equation  and  obtain  the  wing's  circulation 
coefficient  yj  distribution.  With  the  circulation  coefficient 


distribution,  we  can  obtain  the  circulation  coefficient  rate 
of  change  of  the  circulation  coefficient  along  the  wingspan. 

At  the  same  time,  we  can  obtain  the  (4  C^) j  from  the  Yj  and 
after  summation  we  obtain  the  wing 1 s  Cy  lift  coefficient. 

2.  Calculation  of  the  Wing's  Induced  Drag  Coefficient 

If  the  wing  is  in  V  .  then  the  definite  load  distribution 

on  the  wing  can  be  produced.  This  type  of  flow  field  is  called 

the  forward  flow  field.  When  similar  plane  form  wings  are  further 

located  in  -V  ,  if  it  is  able  to  produce  a  similar  load  dis- 
oo 

tribution  to  the  normal  flow  field  then  it  is  called  the  reverse 
flow  field  corresponding  to  the  above  mentioned  forward  flow 
field.  From  reference  [3]  we  know  that  the  wing's  induced  drag 
of  the  forward  flow  field  and  reverse  flow  field  are  equal. 


Fig.  1.  Lattice  Division  and  Location  of  Control  Point 


(1)  JE*» 


(3)  «*«« 


Fig.  2. 


Fig.  2.  The  Forward  Flow  Field,  Reverse  Flow  Field  and  Combined 
Flow  Field 

Key:  1.  Forward  Flow  Field 

2.  Reverse  Flow  Field 

3.  Combined  Flow  Field 

The  combined  flow  field  is  the  superposition  of  the  forward 
flow  field  and  reverse  flow  field.  In  the  combined  flow  field, 
because  the  directions  of  the  forward  flow  field  and  reverse  flow 
field  appended  horse  shoe  vortex  in  the  lattice  are  opposite, 
they  counteract  each  other.  Only  the  two  free  vortex  extending 
into  infinity  are  left  in  front  and  in  back.  At  this  time,  the 
calculation  of  the  wing's  downwash  speed  only  requires  the  cal¬ 
culation  of  the  produced  effects  of  the  free  vortex.  The  in¬ 
duced  drag  of  the  wing  in  the  combined  flow  field  is  twice  that 
when  it  is  either  in  forward  flow  or  reverse  flow.  Therefore, 
the  calculated  induced  drag  must  be  decreased  by  half. 

In  the  combined  flow  field,  the  circulation  coefficient  of 
each  lattice  free  vortex  is  a  known  value  which  is  the  yj  in 
the  advancing  forward  flow  field.  With  the  circulation  coef¬ 
ficient  distribution,  we  can  obtain  the  normal  wash  Dj  of  the 
entire  free  vortex  in  each  control  point  location.  After  the 
summation  of  each  lattice,  we  obtained  the  induced  drag  coef¬ 
ficient  of  the  wing: 

(aC,),D,(AS),  13  > 

i  - 1 

{A  C  ) j  is  the  load  coefficient  in  the  positive  flow  field  and 
P 

is  a  known  value. 

If  the  relative  value  of  A  is  taken  to  indicate  the  Cxj , 
the  induced  drag  form  factor  is: 


3.  Analysis  of  the  Calculation  Method 
In  order  to  test  the  applicability  of  the  calculation 
method,  we  calculated  the  A  values  of  four  types  of  rectangular 
wings  with  different  aspect  ratios  and  compared  them  with  the 
values  of  reference  [  4]  .  The  comparative  values  are  listed  in 
table  1 . 


i  \  ■  2 

X  -  4  J 

X  -  8 

X  >  10 

(1)A<*3C> 

1.000*8 

1.00623 

1 

1.01642 

1.04218 

(2)A(?:itC4n 

1.0007 

i.ooro 

1.0160  | 

1.0120 

Table  1 

Key:  1.  This  paper 

2.  Reference  [4] 

It  can  be  seen  from  the  comparison  in  table  1  that  the 
relative  accuracy  of  the  theoretical  values  found  in  reference 
[4]  are  extremely  close  to  the  calculation  results  of  this 
method  and  that  the  accuracy  is  effective  and  satisfactory. 
Furthermore,  the  use  of  a  constant  roll-angle  to  divide  the 
spanwise  lattice  and  a  relatively  dense  wing  tip  lattice  can 
better  reveal  the  circulation  change  of  the  wing  tip  area.  There 
is  not  a  great  change  in  the  wing  tip  form  and  the  reaction  of 
the  A  value  is  relatively  acute.  Therefore,  when  making  a 
comparison  we  can  use  this  method  to  calculate  the  lift  and 
drag  characteristics  of  wings  with  different  winglets. 
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III.  Discussion  of  the  Calculation  Results 

In  order  to  economize  on  calculation  time,  we  calculated 
50  lattices  on  half  wings  and  M  X  N  was  5x10.  On  709  aircraft, 
each  calculation  of  the  conditions  of  a  winglet  took  about  5-6 
minutes. 

At  first,  we  took  the  measured  wings  with  winglets  shown  in 
fig.  3  as  calculation  examples  and  calculated  the  aspect  ratio 
as  ,7V.  =  6.375.  The  wing  area  used  was  S=2x8x25.5  and  in  an  a=0.1 
(arc  degree)  situation  the  lift  and  drag  characteristics  of 
rectangular  wings  with  different  winglets  were  calculated. 


Fig.  3.  The  Measurements  of  a  Rectangular  Wing  With  a  Winglet 

1 .  The  Lift  and  Drag  Characteristics  of  Rectangular  Wings 
With  Upward  Turning  Winglets 

When  the  upward  turned  angled  on  the  winglet  is  0°,  30°, 
45°,  75°  and  90°,  the  calculated  wing  lift  and  drag  character¬ 
istics  are  as  listed  in  table  2. 


& 

c„ 

!  c.,  i 

i 

A 

0* 

0  129548 

0.009385 

1.018665 

SO* 

0.411652 

0.009117 

1.077522 

45* 

9.386130 

0.008574 

1.151774 

75* 

0.328177 

0.007243 

1.346800 

»«* 

0.329244 

0.007601 

1.401237 

Table  2 
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Table  2  The  Various  $  Values  of  the  Lift  and  Drag  Character¬ 
istic  of  Rectangular  Wings 

It  can  be  seen  from  table  2  that  the  larger  the  S  the  larger  the 
A  value.  A  represents  the  ratio  value  of  the  induced  drag  and 
the  minimum  induced  drag.  The  smaller  this  value  the  better  the 
induced  drag  characteristic  and  the  closer  to  the  minimum  in¬ 
duced  drag.  Therefore,  the  comparison  of  the  S  0°  rectangular 
wing  and  & =0°  condition  not  only  causes  the  block  type  winglet 
upward  turning  to  be  able  to  decrease  the  actual  span  length  of 
the  rectangular  wing  but  the  induced  drag  characteristic  can  also 
have  some  advantages.  The  greater  the  upward  turning  angle  the 
greater  the  difference  of  the  induced  drag  characteristic. 

Ig  order  to  analyze  the  reasons,  we  drew  fig.  4  with  the 
Y  and  dy  of  two  rectangular  wings  being  8  =0°  and  §  =75°  for 


Fig.  4. 


The  y  and  Distribution  of  Rectangular  Wings  of 
8  -0°  and  &  -75* 
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It  can  be  seen  from  fig.  4  that  the  y  distribution  of 
8  =0°  is  relatively  close  to  elliptical  distribution  and  the 
y value  of  the  y distribution  of  {?  =75°in  the  vicinity  of  the 
of  the  wing  and  winglet  boundary  line  has  a  sudden  drop.  More¬ 
over,  the  y  values  in  each  location  were  all  smaller  than  the 
value  when  8  =0°.  Therefore,  the  Cy  value  of  the  8  =75°  was 
smaller  than  when  the  value  was  8  =0°.  Because  the  y  distribu¬ 
tions  are  different,  their  ^  value  the  stronger  the  exerted 
vortex.  When  8  =0°,  there  is  a  strong  exerted  vortex  on  the 
wing  tip  and  when  8  =75°  there  is  not  only  a  strong  exerted 
vortex  on  the  wing  tip  of  the  winglet  but  there  is  also  a 
swelling  of  the  in  the  vicinity  of  the  winglet 's  and  wing's 
boundary  line.  There  is  also  a  relatively  strong  exerted  vortex 
at  this  location.  This  then  increases  the  induced  downwash  on 
the  wing  and  furthermore  causes  the  spoilage  of  the  induced 
drag  characteristics. 

2.  The  Lift  and  Drag  of  Rectangular  wings  With  Patterned 
Winglets 

As  regards  the  lift  and  drag  characteristics  of  rectangular 

wings  with  upward  turned  winglets  discussed  above,  the  winglet 

is  a  whole  flat  plate.  If  we  equally  divide  the  winglet  into 

five  horizontal  narrow  lines  chordwise,  when  each  line  has  a 

different  upward  turned  angle,  then  a  type  of  patterned  winglet 

is  formed.  Starting  calculation  from  the  leading  edge  the  upward 

turned  angle  of  each  narrow  line  on  the  winglet  is^^, 

<*?..,  8. ,  and  The  results  of  the  calculations  of  the  lift 

3  4  5 

and  drag  characteristics  of  rectangular  wings  with  various 
patterned  winglets  are  listed  in  table  3. 
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>4 

(Si.  St.  St.  St.  St) 

c,  i 

c., 

A 

<D 

(0*.0\0*.0*,0\) 

0.429548 

0.009385 

l.OlHujj 

(-50*.  -50 *,70 ',70*. TO*) 

0.453810 

0.012075 

1.174282 

(70*. 70*. 70*.  -50*.  -50*) 

0447306 

0.011373 

1.138361 

(0*.  20*.  40*. 00*. 80*) 

0  521052 

0.015679 

1.156575 

(00*. 60*. 40*.  20*. 0*) 

0.429515 

0.00S649 

0.938963 

(60*.  46*.  30*.  15*.  0*) 

0.445285 

0.009315 

0.940880 

©  ; 

(0‘,  15*.  30*.  45*.  0*) 

0-491092 

0.011530 

0.957503 

IKS 

(0*.  45*.  30*.  15*.  0*) 

0455304 

0.009211 

0.889895 

M 

(0*.  80*.  40°.  20*.  0*) 

0.452066 

0.00913S 

0.895231 

(0*.  75*.  50*.  25*.  0*) 

0.448205 

0.01SS83 

1.552973 

© 

(0*.  30*.  20*.  10*.  0*) 

0.453312 

0.009279 

0.904366 

© 

(0*.  45*.  35*.  25*.  15*) 

0.445599 

0.008798 

0.887387 

(0*.  30*.  0*.  30*.  0*) 

0.455909 

0.000768 

0.941150 

« 

(30 ",  15*.  0*,30*.  15*) 

0.449585 

0.009848 

0.975789 

© 

(20*.  0*.  20",  0*.  20*) 

0.449668 

0.009740 

0.964855 

Table  3  The  Lift  and  Drag  Characteristics  of  Rectangular  Wings 
with  Various  Patterned  Winglets 


From  table  3  we  can  conclude  the  following  few  points: 

(1)  Plan@  is  a  rectangular  wing  of  6*  =0°  and  we  took  it  as 
the  criterion  state. for  comparison.  If  the  wing  tip  section  is 
divided  into  a  patterned  wing  tip,  although  this  can  fundament¬ 
ally  raise  the  value,  yet  some  of  the  <4  values  have  large 
changes  while  some  change  only  slightly. 


(2)  Yet  there  is  the  possibility  of  the  patterned  winglet 
causing  A <  I  for  the  rectangular  wing  and  the  winglet  can  also 
cause  its  induced  drag  to  be  smaller  than  the  minimum  induced 
drag  value  of  the  rectangular  wing  with  an  aspect  ratio  of 
A  =6.375. 


Therefore,  the  rectangular  wings  with  patterned  winglets 
are  different  from  rectangular  wings  with  upward  turned  winglets 
We  can  clearly  improve  the  induced  drag  characteristics  of  the 
rectangular  wings  by  manifesting  the  superiority  of  rectangular 
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wings  with  patterned  winglets. 

(3)  By  comparing  the  15  plane  listed  in  table  3  we  can  see 
that  the ^  in  plan®  (0°,  45°,  30°,  15°  and  0°)  and  plan  (0) 

(0°,  45°,  35° ,  25° ,  15°)  are  relatively  small  and  their  C  values 
are  larger  than  those  of  plan®.  Therefore,  as  for  improving 
the  lift  and  drag  characteristics,  they  are  good  plans  worthy  of 
recommendation.  In  order  to  analyze  the  mechanism  of  the 
patterned  winglet  being  able  to  improve  the  lift  and  drag  char¬ 
acteristics  of  the  wing,  each  line  of  the  y  distribution  for 
plans  (I)  and  (§)  are  drawn  in  fig.  5  and  each  line  of  the  dis¬ 
tribution  for  plan®  is  drawn  in  fig.  6. 


Fig.  5.  The  y  Distribution  for  Each  Line  of  the  Rectangular 
Wing  in  Plans  (I)  and 


dl  Y 

Fig.  6.  The  -r-  Distribution  for  Each  Line  of  the  Rectangular 
Wingain  Plan<® 

It  can  be  seen  from  fig.  5  that  each  line  of  the  y  distri¬ 
bution  in  plans  and  (jp  are  noticeably  different.  For  the  y 
values  in  plan  besides  the  winglet  section  having  a  sudden 
drop  smaller  than  plan^).  the  y  values  of  s=2,3,4,5  are  all 
greater  than  plan®.  There  is  relatively  large  swelling  in  the 
winglet  section  and  therefore  the  value  of  plan®  is  greater 
than  the  value  of  plan  Q).  Each  line  of  the  y distribution  in 
plan  CD  is  progressively  decreasing  and  therefore  each  line  only 
has  the  ^  .maximum  value  in  the  wing  tip  area.  However,  there 
is  a  relatively  strong  exerted  vortex  on  the  wing  tip.  The  ^ 

of  plan are  all  negative  values  in  the  four  lines  of  s=2,3,4,5 
and  moreover  the  minimum  value  point  of  the  existing  — f 

shown  at  these  points  to  have  a  relatively  strong  reverse  ex¬ 
erted  vortex.  The  reverse  exerted  vortex  produced  induced  upwash 
on  the  wing  and  thus  causes  a  decrease  of  the  induced  drag. 

It  can  clearly  be  seen  from  fig.  6  that  plan®  still  has  a 
positive  maximum  value  in  each  line  of  the  pointed  end  of  the 
winglet  yet  because  the  s=2,3,4  lines  of  the  winglet  turn  up¬ 
wards  one  degree  these  forward  exerted  vortex  are  not  again  in 
the  x-y  plane  but  shift  a  distance  in  the  z  direction.  This 
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also  decreased  the  z  direction  induced  downwash  weight  on  the 
wing  and  was  partially  favorable  to  decreasing  the  induced  drag. 

In  order  to  further  observe  the  wing  tip  exerted  vortex  of 
plan  we  drew  schematic  fig.  7. 


Fig.  7.  Schematic  Chart  of  the  Wing's  Wing  Tip  Exerted  Vortex 
of  Plan 

It  can  be  seen  from  the  chart  that  in  the  wing 1 s  wing  tip 
and  its  vicinity  there  not  only  appears  a  relatively  strong 
forward  exerted  vortex  but  there  also  appears  a  relatively 
strong  reverse  exerted  vortex.  At  the  same  time,  the  winglet 
has  dispersion  and  upwards  shift  of  the  vortex.  Therefore,  the 
winglet  in  plan  (D  has  noticeable  vortex  elimination  and  is  able 
to  improve  the  induced  drag  characteristics  of  the  wing. 

3.  The  Lift  and  Drag  Characteristic  of  Upward  Turning 
and  Downward  Turning  Winglets 

In  order  to  compare  the  effects  of  the  upward  turning  and 
downward  turning  of  winglets  on  the  lift  and  drag  character¬ 
istics  of  wings  we  made  the  following  few  group  plan  comparative 
calculations  (see  table  4) . 
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»  : 

( 6  i.  ft  i>  3  j.  6  «,  &  s) 

(6  i»  6  i»  &  s.»  &  4«  &  5) 

<30*»  30*,  30",  30*,  30*) 

<  -  30*.  -  30*,  -  30*.  -  SO*,  -  30*) 

© 

<75*,75*,75*.75°,75*> 

(-75*.  -75*.  -75*.  -  75*,  -75*) 

® 

<80*,60‘,  10*,20“,0*) 

<-80*.  -  SO*.  -40*.  -  20*,  0*) 

(0°,  20*,  40**60*,  80") 

<0\  -30*.  -40*,  -60*,  - 80 *> 

® 

<0*,  15*, 30*. 45*,  0*) 

(0*,  -  15*.  -30*,  -  45*,  0*) 

<0*.  -45*.  -  35*.  25*  15*3  1 

,0*  15*  .  35  *.  -  23  ’,  -  1-3*) 

Table  4  Comparative  Plan  of  the  Upward  Turning  and  Downward 
Turning  of  the  Winglet 


Results  showed  that  the  calculated  Cy,  Cxj  and  h  values  of 
the  upward  turning  plan  and  downward  turning  plan  were  totally 
identical-  This  explained  that  only  if  there  was  symmetrical 
upward  turning  then  the  lift  and  drag  characteristics  would  be 
identical  to  the  upward  turning.  It  can  be  seen  from  table  4 
that  the  combined  upward  and  downward  turning  of  the  wing  tip 
only  needed  an  unchanging  turning  angle  increase  or  decrease. 

At  the  same  time,  the  included  angle  adjacent  to  the  winglet 
does  not  change  so  that  its  lift  and  drag  characteristics  also 
do  not  change. 


4.  The  Induced  Drag  Characteristics  of  Different  Sized 
Winglets 

The  previously  calculated  lift  and  drag  characteristic 
values  were  all  concerned  with  the  reference  area  as  the  wing 
area  after  the  winglet  was  flattened  for  the  a=0.1  (arc  degree). 
Now  we  will  discuss  another  situation,  that  is,  using  the 
A.  =4.5  rectangular  wing  of  plan  as  the  basic  wing.  When 
without  a  winglet  the  wing's  reference  area  is  S=2x8xl8.  If 
different  sized  and  different  patterned  winglets  are  added  on 
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to  this  basic  wing  and  use  plans  and  <§,  although  the 

lengths  of  each  of  the  added  winglets  of  these  four  plans  is 
different  yet  their  actual  projections  on  the  half  wingspan 
length  on  the  y  axis  are  all  close  to  20.8  (see  fig.  8). 


Fig.  8.  The  Plan  of  Different  Sized  Winglets 

If  plans  ($).  <^)  and  (§)  all  use  the  reference  area  of  plan  (£,  their 
aspect  ratios  all  use  A  =4.5  and  the  same  values  of  plan  are 
comparable  to  their  C  .  values.  This  then  requires  that  the  orig- 

Aj 

inal  lift  and  drag  values  of  plans(J)/  @  and  be  converted 
into  conversion  values.  The  conversion  values  of  each  plan  are 
listed  in  table  5. 


JU 

<6 1>  & i»  6 «>  6 »> 

<**> 

(2)"' 

C, 

caf 

A 

A  C.,(H> 

C 2) 

0 

!  »•! 

0.370347 

0.010265 

1.008454 

0 

(0**0*.  0*i  0**0")  j 

2  2 

j  0.085214  . 

0.370347 

0  008172 

0.8029t8 

20.3818 

<0\45\30MS\0*> 

2.2 

|  0.081428  { 

0.379347 

0.007865 

0.77273 

23.3719 

(SS.S’.SO.S’.SO.S’.SS.S’.SS.S*) 

4.2 

j  0.081335  | 

0.379347 

0.007675 

0.75404 

25.2282 

7.5 

|  0.083832  i 

0.379347 

0.006755 

0.663623 

34.1941 

Table  5  The  Conversion  Values  of  Each  Plan 

Key:  l.  winglet  length  (centimeters) 
2.  a  (arc  degree) 
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It  can  be  seen  from  the  comparison  in  table  5  that: 


(1)  The  added  2.2  length  of  patterned  winglet  plan^)has  a 
23.37%  decrease  in  induced  drag  as  compared  to  plan which  does 
not  add  a  winglet.  This  is  better  than  the  aspect  ratio  of 

plan Q which  has  an  induced  drag  decrease  of  20.38%. 

(2)  The  7.5  added  length  simple  upward  turning  winglet  plan(f 
reduced  the  induced  drag  by  34.19%  more  than  plan  (p and  is  better 
than  plans  (J),  ©and(^l  Although  these  plans  have  the  most  de¬ 
creases  in  induced  drag  yet  the  winglet s'  measurements  are  large, 
are  heavy  in  weight  and  the  winglets’  y  direction  lateral  force 
can  greater  enlarge  the  bending  moment  of  the  wing  root.  When 
selecting  plans,  we  should  synthesize  each  aspect ' s  good  and  bad 
points  and  make  the  choice  by  means  of  an  overall  comparison. 
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THE  CALCULATION  OF  LIFT  AND  DRAG 
CHARACTERISTICS  OF  SUBSONIC  WINGS  WITH 

WINGLETS 

Zhou  Renliang 

( Nanjing  Aeronautical  Institute') 

Abstract 

In  order  to  calculate  the  lifts  of  rectangular  wings  with  different 
winglets  at  subsonic  speeds,  we  have  adopted  the  finite  element  method 
which  divides  spanwise  lattice  and  determines  spanwise  locations  of  control 
points  by  means  of  a  constant  roll-angle  method.  The  induced  drags  are  also 
calculated  by  using  combined  flow  field  method. 

As  the  results  of  calculating  various  configuration  of  a  winglet,  we  have 
found  out  some  rules  affecting  the  lift  and  drag  characteristics  of  wings  with 
winglets  and  picked  out  a  favourable  configuration  from  them.  The  aerodyna¬ 
mic  mechanism  of  winglets  is  also  discussed. 


A  NEW  MODEL  FOR  PREDICTING  OVERLOAD  RETARDATION  EFFECT  IN 
FATIGUE  CRACK  PROPAGATION 

by  He  Qingzhi 

(Beijing  Institute  of  Aeronautics  and  Astronautics) 

Abstract 

Based  on  the  concept  of  effective  residual  compressive 
stress  and  considering  the  residual  compressive  stress  relax¬ 
ation  in  basic  circulation,  we  calculated  a  formula  for  the 
propagation  speed  of  cracks  after  overload: 


fda\ 

...  /  da  \ 

UW, 

i 

M 

I 

(da) 

In  the  formula, (dn)  is  the  propagation  speed  of  cracks  without 

kap 

the  effects  of  overload;  u  =  1  +  (1-a)  A.  -  (1-a)  x  4Kb;  the 
stress  relaxation  coefficient  is: 


a 


1  ~  MCt 

r.-l 


The  method  in  this  paper  was  used  to  calculate  some  specimens 
of  LY12 ,  2024,  7075  aluminum  alloys,  Ti-6A1-4V  titanium  alloy  and 
a  stiffened  plate  under  different  loading  conditions.  By  putting 
in  the  coefficient  n=4 ,  we  discovered  that  the  calculated  re¬ 
sults  were  in  good  agreement  with  the  test  results. 


52 


I.  Preface 

Given  an  alternating  circulation  loading  component,  after 
exerting  single  or  multiple  high  peak  loads  (overload) ,  the 
crack  propagation  speed  decreases  noticeably,  after  a  period 
of  circulation,  the  crack  propagation  speed  then  gradually 
returns  to  normal.  Fig.  1  (a)  shows  loading  conditions  and 
fig .  1  (b)  is  a  schematic  chart  showing  the  crack  propagation 

characteristic  points  after  overloading.  The  above  mentioned 
phenomenon  was  previously  given  attention  and  a  great 
deal  of  research  work  was  carried  out. [1-6]  Many  theories  were 
proposed  to  clarify  the  physical  quality  of  this  phenomenon. 
Among  them,  the  most  notable  are  the  theory  of  the  residual  com¬ 
pressive  stress  in  the  plastic  zone[l,3!and  the  crack  closure 
theory  [  9] .  Wheelerf  7  ]and  Wlllenborg  [8]  separately  proposed  models 
for  calculating  the  retardation  period.  These  models  which  were 
used  to  estimate  the  load  spectrum  of  the  retardation  period 
under  loading  were  relatively  successful,  yet  they  did  not 
predict  the  retardation  effects  of  single  overload  very  satis¬ 
factorily. 

This  paper  uses  the  concept  of  effective  residual  compres¬ 
sive  stress  to  propose  a  new  model  which  can  accurately  predict 
the  retardation  effects  after  a  single  elongated  overload. 


Fig.  1. 
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Key:  (a)  Loading  conditions  of  overload 

(b)  Schematic  chart  of  the  crack  propagation 
characteristic  points  after  overloading 

II.  Model  and  Calculation  Formula 

After  a  single  elongated  overload,  there  is  residual  com¬ 
pressive  stress  in  the  plastic  zone  produced  by  the  overload. 
Based  on  the  hypothesis  found  in  reference  [8],  this  compres¬ 
sive  stress  is  equivalent  to  an  externally  added  compressive 
stress  (f  on  the  specimen: 

(  i  ) 


In  the  formula,  0^  is  the  circulation  stress  peak  value  required 
for  this  instantaneously  caused  retardation  effect  loss.  Its 
value  is: 


/Cn  y/R^-Aa 
Y(a) 


(2) 


See  the  appendix  for  the  deduction  of  formula  (2)  and  for  the  de¬ 
finitions  of  each  of  the  symbols. 


Taking  the  hypothesis  a  step  further,  in  basic  circulation 
(low  load  circulation)  when  the  stress  drops  from  the  peak  value 
<y]binax  to  the  valley  Cf^^,  there  are  changes  in  the  stress  field 
of  the  crack  tips.  At  the  same  time,  the  crack  propagates  forward 
a  micro- 8  a  and  there  is  local  relaxation  of  the  residual  com¬ 
pressive  stress  in  the  overload  plastic  zone;  because  of  this, 
the  equivalent  externally  added  compressive  stress  also  decreases. 


that  is,  it  drops  from  (7  to  ao*  .  It  is 

r  r 

(or).,.-a(o.rO  (3) 

. . 

In  the  formula,  a  is  the  compressive  stress  relaxation  coef¬ 
ficient  whose  value  will  be  determined  below. 

Because  of  the  action  of  an  equivalent  residual  compressive 
stress,  the  peak  value  load  and  valley  load  effective  values  of 
the  basic  circulation  (see  chart  1(a))  are  separately: 

(Ofcn .  )  ,//  — °»m.x  ~°r 

(4) 


We  can  obtain  the  effective-  stress  amplitude  value  from  formulas 
(3)  and  (4)  : 


(A  — 


(o*mi „)«//  =  A<*»-  (  1  -  °  )  (°*f  a‘-* * 


) 


(5) 


In  the  formula,  A  0L  -  Cf  -  O',  .  Formula  (5)  can  also  be 

d  bmax  bmm . 


written  as: 


(A^»)«»"AJCi-  (1  — ®)  (^«»_^(«u) 


(6) 


In  the  formula:  (aAT»).,/- k( a  )(Ao,).„ 

AAT»*K(«)Aa, 

K(  a  )°»«« 

Y  (a)  is  the  coefficient  of  the  crack's  characteristic  length 
a  and  the  geometrical  shape  of  the  component. 
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When  the  two  ends  of  formula  (6)  are  divided  by  4k  ,  we 

b 

obtain : 


==(A^£ZLatl+(1_a)x_(i  — a)  Kg_ 


(7) 


aa:» 

In  the  formula: 

x  =  ~-=  (1) 

A/C*  1  /?  (and)  a»-«  A*-» 

1.  is  the  stress  ratio 

It  is  very  easy  to  prove  that  when  1,  u  is  the  reduced 
coefficient  of  the  stress  strength  factor  amplitude  value  after 
overload . 

Formula  (7)  can  be  written  as: 


(aA'6)„,  =  u\/v 


C «  ) 


If  the  crack  propagation  rate  is  expressed  by  the  Paris  formula, 
then 

S8-  ci  (^K)n 

It  can  be  seen  that  the  effective  value  u(/J 

of  the  crack  propagation's  stress  strength  factor  amplitude  is 
caused  after  overload  and  therefore  the  cracking  propagation 
speed  will  decrease  after  overload. 

Below  we  will  seek  the  compressive  stress  relaxation  coef¬ 
ficient  a. 


It  can  be  seen  from  formula  (22)  that  when  &  a  =  0  (.that  is, 

just  after  overload) ,  K  has  a  maximum  value. 

ap 


(*.,) 


(9) 
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We  know  from  formula  (6)  that  at  this  time  (4  K^)e£f  has  a 
minimum  value 


-AiC»Ci  -(1  -«)(••-  D3 


(10) 


In  the  formula. 


AC, 
r  '44 


..  .  ,  F.N.l 

=  -  is  the  overload  ratio 

4)Ai> 


Many  tests  have  proved  that  when  overload  ratio  r  is  greater 

than  a  certain  critical  value  r  (the  r  value  is  related  to  the 

c 

data  and  stress  ratio  R;  further,  it  increases  in  accordance  with 
the  increases  of  R) ,  crack  propagation  will  be  completely  re- 
tarded[17,18](F'N*1)  .  At  this  time. 


(AK Ai<*C  1  -  <  1  -  1  )) 


F.N.l:  In  references  [17]  and  [18]  the  overload  ratio  is 

y  K  ^ 

defined  as  r  *  lmax  .  Naturally,  between  r  and  r' 

Dmax  r'-R 

there  is  the  following  relationship,  r  =  — ~R. 
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Naturally,  if 


(  A  A  *)  A  A 


(11) 


crack  propagation  will  be  completely  retarded;  A K  ^  is  the 
"threshold  value"  of  the  fatigue  crack  propagation.  Using  the 
conditions  expressed  in  formula  (11) ,  we  can  obtain: 


a 


A  K,t 
~  A  A, 

r,~  1 


(12) 


When  formula  (12)  is  substituted  into  formulas  (7)  and  (8) ,  we 
obtain : 


(A  AO 


•tr 


.  AA,»  i  AAj* 

1  aA7~  '  1  A  A* 

r.-  l  *  r.-  1 


AA»=h(AA*) 


(13) 


During  the  retardation  period  following  overload,  the  crack  pro¬ 
pagation  speed  is: 

“0  (2)jjj  (AA0,«<AA,. 

®(3)  (*^)f  -«*  *  (14) 

Key:  1.  If 

2.  If 

3.  Or 
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In  the  formula, 


is  the  crack  propagation  speed  corres¬ 
ponding  to  &  Kb  without  being  affected  by  overload  and  n  is  a 
coef ficient . 

Integral  formula  (14)  can  obtain  the  circulation  number 
Nj  required  to  penetrate  the  overload  retardation  zone. 


.V 


d&a 


(15) 


ac  is  the  length  of  the  overload  retardation  zone  and  its 
value  can  be  obtained  by  putting  in  u  =  1  from  formula  (7) : 


(16) 


The  comparison  of  retardation  zone  length  A  ac  and  crack 

characteristic  length  a  is  very  small  and  when  A  CT^  is  the 

constant,  we  can  consider  the  A  K,  and  as  the  constants 

b  (dN) o 

At  this  time,  formula  (15)  can  be  written  as: 


Below  we  will  study  the  changes  of  the  crack  speed  after 
overload.  It  can  be  seen  from  formulas  (7)  and  (22)  that  the  u 
value  increases  in  accordance  with  the  increases  of  the  Zl  a  and 


when  it  reaches 
(da) 

(dN) 


Z)  a  =  A 


u  =  1.  The  minimum  values  of  u  and 


are : 


l  -(1  -<*)(*■-  1) 


(18) 


and 


(19) 


III.  Results  and  Discussion  of  Calculations 


Applying  the  previously  deduced  formulas,  we  separately 
calculated  the  overload  retardation  period  of  some  specimens  of 
LY12 ,  2024-T3,  2024-T351,  7075-T6,  7075-T73  aluminum  alloys, 
Ti-6Al-AV  titanium  alloy  and  a  stiffened  plate  under  different 
loading  conditions.  The  Zl  and  rc  data  are  separately  taken 
from  references  (12,13,18,19,20). 


In  the  calculations,  we  used  the  coefficient  n  =  4.  The 

calculation  results  of  N*  are  listed  in  table  1  and  the  calcula- 
(da)  Jn  da) 


tion  results  of 


(dN) 


(<Hi  n  •  -o 

tables  list  the  corresponding  test  values. 


are  listed  in  table  2.  Both 
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w  # 

(1) 

A/C» 

MN 

— 

R 

««& 

(3) 

r 

«*/* 

c®(7 

.  K‘m  (5) 

TFJ - 

A  Mil  ft) 

N* (S-.fi)  (9) 

(»*«*> 

3. 
m  * 

lBk  ft  a 

LYl2 

9.4 

0.1 

2.0 

16.5xl0-»“'> 

1.0 

30000 

S2000r333 
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<4> 
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0.1 

2.0 
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vO 
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1.5 
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16.5 
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id. 8 

0 

1.5 
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14 

0 
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150001133 
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! 

14 

0 
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7.6 x  10*®f3®3 
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0.8 

14 

0 
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SOOOOC*33 
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Ti-6A1-4V 
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0.1 

1.78 

17. 8x  10*®13®3 

1.0 
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17.4 
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1 

i 
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1 

1.3 

7075-T73  ! 

7.81 
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2.11 
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250001 143 
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<1.6> 
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i 
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Table  1 
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Table  t 


Key:  1.  Material  (thickness  in  millimeters) 

2.  A  model  stiffened  plate 

3.  Overload  ratio 

4.  Millimeters/cycle 

5.  (Cycles) 

6.  (Calculation) 

7.  Calculated  value 

8.  Test  value 
9  *  Test 

10.  (J)  Calculation  based  on  data  from  reference  [10] 

11.  @  Data  of  503  tests  at  the  Beijing  Institute  of 
Aeronautics  and  Astronautics 


tf  ft 

(1) 

0*K**> 

Ak* 

R 

(2) 

r 

(zr)„h./(d7r)0 
(3)  <«•#> 

( dN  )mi»  / (dN  )„ 

(4)  (*»)© 

(W-L* 

(awL. 

13.1 

0.01 

D 

1.65% 

1.55 

14. » 

0.01 

1.59% 

1.38 

2024-T351 

18.0 

0.01 

tMH 

1.45% 

1.26 

(10) 

13.5 

0.1 

■n 

3.45% 

1.5  ?» 

1.55 

15.7 

0.1 

Bn 

2.09% 

1.32% 

1.58 

20.2 

0.1 

Bn 

1.66% 

2.5*% 

0.66 

2024-T3 

(3.2) 

19.0 

_ 

0.05 

D 

1.67% 

1.5% 

1.11 

(7)®  *m*»**XKa*&*l. 


Table  2 

Key: 


1.  Material  (thickness  in  millimeters) 

2.  Overload  ratio 

3.  (Calculation) 

4.  (Test) 

5.  (Calculation) 

6.  jTest) 

7.  Q  Data  taken  from  references  [14]  and  [4] 


The  comparison  of  N*  (calculation)  and  N*  (test)  is  shown 
in  fig.  2. 


Fig.  2. 

Key: 


1.  (Calculation) 

2.  (Test) 

3.  Stiffened  plate 


After  the  above  analysis  and  calculations,  we  can  make  the 
following  conclusions: 

(1)  Coefficient  n  =  4  can  be  used  for  the  LY12,  2024,  7075 
aluminum  alloys  and  Ti-6Al-4V  titanium  alloy; 


(2)  The  overload  retardation  period  N*  obtained  from  the 


calculations  of  the  proposed  model  used  in  this  paper  were 

in  good  agreement  with  the  test  results.  When  considering  the 

dispersity  of  the  crack  propagation  speed  and  the  N* 

1  o 

test  data,  we  found  that  using  the  method  in  this  paper  to  cal¬ 
culate  the  retardation  period  was  quite  satisfactory? 

(3)  The  crack  minimum  propagation  speed  after  over- 

'  min 

load  calculated  by  means  of  this  method  was  generally  slightly 
greater  than  the  test  values.  Moreover,  when  occuring  just  after 
overload,  this  speed  was  not  in  agreement  with  the  test  results. 

Appendix 

The  deduction  of  formula  (2) : 

Based  on  the  concept  proposed  in  reference  [8],  we  have 
(see  fig.  3) : 


<»i  +  #>.,s=ao  +  AJ>1 


(20) 


Fig.  3. 
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In  the  formula,  R 


yi 


1 

C7T 


Cl 


«cf 


max 


ys 


J2 

) 


Plastic  zone  diameter  pro¬ 
duced  by  overload; 


R  _  _i_  (Kap)  2 
yap-  Clf  (tfys) 


-  the  plastic  zone  diameter 
derived  from  the  added  cir¬ 
culation  stress  peak  value 

<T _ of  a  certain  instantan- 

ap 

eous  and  hypothetically 
caused  overload  effect  loss 


C  -  the  constant  of  the  Irwin 
plastic  zone 


C  =  1  is  the  situation  for  plane  stress; 

C  =  3  is  the  situation  for  plane  strain; 

C  =  3  is  the  situation  for  the  combined  model,  s  is 

l+2s 

the  proportion  occupied  by  the  inclined  fracture  on 
the  break  fracture  (shearing  model) . 


From  formula  (20) : 
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is  the  crack  propagation  quantity  after  over- 


If  A  a  =  a.  -  a 
i  o 

load ,  the 


Out/  Cn y/  Rti  — 


(21) 


or 


„  o^y/Cn  »//?,,  —  Aa 

°"  K(a) - ujf - 


(22) 
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A  NEW  MODEL  FOR  PREDICTING  OVERLOAD 
RETARDATION  EFFECT  IN  FATIGUE 
CRACK  PROPAGATION 


He  Qingzhi 

( Beijing  Institute  of  Aeronautics  and  Astronautics ) 

Abstract 


The  effect  of  the  residual  compressive  stress  in  the  overload  plastic  zone 
is  equivalent  to  an  effective  stress  Or,  and  as  proposed  by  Willenborg 

where 


j/ H x  Au 

Y(a) 


In  the  present  paper  it  is  further  assumed  that  or  ;s  not  constant  during 
one  baseline  cycle,  but  it  varies  from  o,  at  <*»  mtM  to  «ar  at  o*  mim,  «  is  a  rela¬ 
xation  factor.  Then 

••»)»//*  Ot  o, 

(®t  mi »)»//  3  Ok  min 

so  the  effective  stress  range  (A <*»)./»  during  the  baseline  cycle  is 
(A°i)««*  (ok  (cr»  Act*—  (  l  —  o  )  (a.,— a* 

Or  alternatively  the  above  equation  may  be  written  in  the  following  form 
(AfC*)*//**  (  1  —  o  —  KimM) 

It  is  easy  to  show  that  the  maximum  value  of  Kmp  during  the  transient 
period  following  the  overload  will  be 
i  m*m  at  Aa  =  0,  so 

(aa:»),m«,.“Aa:,c  i  -( 1  -  <*)(  r  -  i  )j 
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where  the  overload  ratio  r  = 

Many  investigators  have  found  that  when  the  overload  rate  increases  to 
some  critical  value  r„  the  crack  will  cease  to  propagate  at  all  after  the  over¬ 
load.  It  is  reasonable  to  assume  that  when  r  reaches  r„  (A u.  will  de¬ 
crease  to  a  value  „u=  aA,*  so  the  crack  will  cease  to  propagate. 

Using  this  condition,  we  get 

AJC»Cl  -(1  -a)(r.-  D3-AA',* 


or 


a 


A  Ktk 

Z_  WL 

r.—  1 


If  the  crack  propagation  rate  is  expressed  by  Paris  formula 

^-C,<4  KiT 

then  the  reduced  rate  of  crack  propagation  following  an  overload  will  be 

where  «•  l+(l— <*)X— (l  — 

1  JSH. 

K  &K„ 

{4k\  - c '^K‘r  the  unretarded  rate  of  crack  propagation  corresponding 
to  (A*.) . 

The  retarded  number  of  cycle  N*  can  be  calculated  by 


N.,  f* _ dAa 

D  .  .( d<r 


I 


where  A««  is  the  length  o  he  retarded  zone. 

The  retarded  number  of  cycle  of  some  specimens  of  LY12,  2024,7075  alu¬ 
minium  alloys,  Ti-6A!-4V  titanium  alloy  and  a  stiffened  plate  under  diffe¬ 
rent  loading  conditions  are  predicted  by  the  present  method.  By  putting  the 
coefficient  n  *  4 ,  we  find  that  the  agreement  between  the  predicted  values 
and  the  tested  values  is  good,  as  shown  in  fig.  2  of  this  paper. 
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J-INTEGRAL  EXPERIMENTAL  CALIBRATION  OF  SHEET  SPECIMENS  WITH 
SINGLE  EDGE  NOTCH 


Lui  Ligeng,  Chen  Xianxi  and  Cai  Qigong 
(Central  Iron  and  Steel  Research  Institute) 

Zheng  Minzhung 

(Aircraft  Strength  Research  Institute) 


Abstract 

The  J-integral  of  a  sheet  notched  specimen  has  a  very  simple 
form  when  the  notched  surface  has  been  taken  as  the  integral 
contour,  i.e. 


/  =  J  Wdy* 


p/2 


*12 


fT(B)  PcosddO 


It  can  be  further  simplified  into: 

J=APW 

o 


where 


J*/2 

o  «P  (  8  )d& 
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and 


( e }  =  -^p-c°s 9 


For  two  types  of  specimens  with  single  edge  notch  (the 
notch  root  curve  radius  is  divided  into  20  millimeters  and 
9  millimeters) ,  this  paper  determined  the  strain  distribution 
£ (#)  of  the  notched  surface,  then  after  stress-strain  curve 
was  converted  to  form  deformation  work  density  distribution 
W  (0)  there  was  finally  numerical  integration  and  we  obtained 
coefficient  A.  The  test  results  were  compared  with  the  linear 
elasticity  fracture  mechanics  method  of  Weiss's  correction  of 
the  plastic  zone. 


I.  Preface 

When  analyzing  the  strength  and  predicting  the  lifespan  of 
the  notch  specimen,  it  is  necessary  for  people  to  determine  the 
relationship  between  the  outer  load,  specimen  and  notch  geo¬ 
metrical  parameters,  and  the  notch  root's  maximum  strain 
the  maximum  stress-strain  product  £q  <Xq  and  the  maximum  deforma¬ 
tion  work  density[L-3]W0=  ^ <fd£.  It  is  well  known  that  when  the 
notch  root  is  located  in  an  elastic  sphere,  £  =K.  £  , 

2  A2  2J2  9 

£  a  =k:  a  and  W  £  <f  =i5  Kf  _o  •  In  this  equation,  6  and 
°  °  ^  p  o  *  o  o  trr*  n  g 


cf  are  the  distant  mean  strain  and  stress,  and  K  is  the  theor- 

etically  elastic  concentration  coefficient.  After  the  notch  root 

entered  into  an  elastoplastic  state,  it  went  through  the  commonly 

2 

used  Neuber  relationship,  that  is  K.K.  =K  ,  and  obtained 

o  =  o  t 


*o3o 


-  £t<V  K* 


Based  on  the  J-integral  full  quantity  theory,  in  reference 
[4]  Cai  Qigong  proved  that,  under  the  conditions  of  a  hardened 
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material  of  an  assumed  power  and  when  the  notch  root  strain 

distribution  satisfies  the  £(8)  =  £  cosm9,  the  Neuber  rela- 

o 

tionship  is  still  established  under  a  type  I  load.  He  also  cal¬ 
culated  the  A  in  J/ (^AWq  when  in  the  0.8-1. 2  range. 

Therefore,  after  J-integral  direct  experimental  calibra¬ 
tion  of  the  notch  specimen,  we  measured  the  £(S)  and  W(8)  and 
experimentally  determined  coefficient  A.  Thus,  we  used  the 
fracture  mechanics  parameters  to  describe  the  maximum  deforma¬ 
tion  work  density  W  of  the  notch  root. 


II.  Specimen  and  Experimental  Method 
1 .  Specimen 

The  specimen  was  taken  from  hot  rolled  sheet  of  aluminum 
alloy  LY12  with  a  50  millimeter  thickness  of  Cfg  2=H*2  kilo- 
grams/millimeter^ .  The  measurements  of  the  sheet  speciment  with 
a  single  edge  notch  was:  length  400  millimeters,  width  150 
millimeters,  thickness  30  millimeters.  See  fig.  1  for  the  speci¬ 
men  number  and  notch  measurements. 


<T  n 
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?ig.  1.  Measurement  of  Specimen  Notch 
KeY:  1.  Specimen  number 
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r 


The  following  formula  can  be  used  to  indicate  the  stress- 
strain  curve  of  the  measured  LY-12. 


t 


=  — — +  ao' 
E 


(  1  ) 


In  the  formula,  a=3. 537x10  n=9.25  and  E=7000  kilograms/ 

millimeter^ . 

2 .  Test  Method 

The  test  was  carried  out  on  a  200  ton  East  German  ZDM200TPU 
type  material  teit  machine.  The  two  ends  of  the  specimen  were 
each  fined  on  a  plate  with  6  bolts  and  the  plate  used  dowels 
to  connect  the  connecting  rod  to  the  test  machine.  In  the  test, 
specimen  C^gb  was  directly  clamped. 

There  were  lxl  millimeter  and  2x3  millimeter  resistance 
strain  pieces  pasted  on  the  inner  surface  of  the  specimen  notch. 
Their  model  numbers  and  resistance  values  were  PBJ^-A.^  (R=60  ohm) 
and  PBJ^-A^  (R=120  ohm) .  The  resistance  strain  pieces  were  sep¬ 
arately  pasted  on  the  0°,  10°,  20°,  30°,  45°,  60°,  75°  and  90° 
curvature  central  arguments  of  the  corresponding  notches  (see 
fig.  2.) 


71 


Jfi  i  t  ^2Qb 


Fig.  2.  Developed  Distribution  Chart  of  the  Strain  Piece  on 
the  Inner  Surface  of  the  Notch 

Key:  1.  Location  of  pasted  strain  piece 

At  the  time  of  the  test,  there  was  gradual  loading  so  as 
to  maintain  the  load  and  measure  the  strain  value  up  until  the 
strain  piece  lost  effectiveness  when  it  reached  the  maximum 
strain  area  at  about  9000  microstrain. 

III.  Test  Results 

1.  The  Strain  Distribution  of  the  Notch  Root 

Figs.  3-5  present  the  distribution  of  the  strain  along  the 
notch  surface  under  different  load  levels.  In  two  specimens 
with  relatively  large  curvature  radii,  there  was  a  strain  dis¬ 
tribution  form  similar  to:  a  specimen  with  a  curvature  radius 
of  p=9  millimeters  had  an  obvious  widening  of  the  strain  concen¬ 
tration  region. 


Fig.  3 


.  The  Strain  Distribution  of  Specimen  Q  Under  Different 
Load  Levels  a 

Key:  1.  Microstrain  (6) 

2.  Angle  0  (degrees 

3.  Load  (tons) 


73 


Fig.  4 


Fig.  5 


.  The  Strain  Distribution  of  Specimen  Q  „  Under  Different 
Load  Levels 

Key:  1.  Microstrain  (9) 

2.  Angle  6  (degrees) 

3.  Load  (tons) 
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Fig.  5.  The  Strain  Distribution  of  Specimen  Q  Under  Different 
Load  Levels  ^ 

Key:  1.  Microstrain  (©) 

2.  Angle  8  (degrees) 

3.  Load  (tons) 


The  tests  also  found  that  when  three  specimens  were  in  the 
0=90°  area  and  the  A, 8  and  C  areas  (see  Fig.  1.)  on  the  outer 
side  of  the  notch,  they  all  appeared  in  the  compressive  stress 
region  and  the  strains  were  all  negative  values. 

2.  The  Deformation  Work  Density  Distribution  of 
the  Notch  Root 

The  conditions  of  the  plane  stress  was  (j=d=0 .  Because  of 

r  z 

this,  the  deformation  work  density  was: 


r  t  (8> 

*r(0)  =  J  o  adt 


( ::  > 


In  the  formula  and  afterwards,  the  o  and  t  were  separately  in¬ 
dicated  as  cf  and  £  . 

e  • 

d 

The  conditions  of  the  plane  strains  were  C  =0,  6  =  °  * 

r  z  2 

€  =0  and  t  =  -Hr  -If  the  material's  equivalent  stress  and  strain 
curves  can  oe  approximately  indicated  by  the  following  formula 


(  3  ) 
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then 


PF(0) 


( 4  ) 


The  deformation  work  density  under  mixed  conditions  can  be 
written  as: 


^(e)=/(m)j''<9>ac/e 


Therefore,  based  on  the  material's  stress-strain  relation¬ 
ship,  we  can  convert  the  strain  distribution  £ (0)  of  the  notch 
root  into  the  equivalent  deformation  work  density  distribution 
which  is  W(9)/f(m)  (see  fig.  6.). 


Fig.  6. 
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Fig.  6.  The  Deformation  Work  Density  Distribution  of  Special 

^20a  Under  Different  Load  Levels 

2 

Key:  1.  Kilograms/millimeter 

2.  Angle  Q  (degrees) 

3.  Load  (tons) 


In  the  specimen  with  a  p=20  millimeter  notch,  the  W(8)  can 
be  approximately  indicated  as : 


(V(e)  =  fVtcoa*Q 


(6) 


For  the  specimen  with  a  p=9  millimeter  notch,  then: 


(7  ) 


3.  Determining  Coefficient  A 

Based  on  the  definition  and  characteristics  of  the  j-inte- 
gral ,  we  can  take  the  notch  surface  as  the  integral  contour  at 
which  time: 


1 


-ir— n 


e  jpcwede-^pH'', 


(x) 
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In  the  formula 


A  =  2  J“/a  4>(  0 


In  the  formula,  ¥  (©)=  ^ ^  cos  8.  When  formulas  (6)  and  (7) 

w  o 

are  separately  substituted  into  formula  (9)  we  can  analytically 
obtain  coefficient  A. 

A=0 . 91  when  p=20  millimeters  or  a/p=2 
A=1 . 18  when  p=9  millimeters  or  a/p=4.4 

Based  on  fig.  6,  we  can  also  convert  the  W(8)/f(m)  figure  into  a 
¥  (9)  figure  (see  fig. 7.).  We  can  obtain  the  A  value  by  num¬ 
erical  integration  (see  table  1) . 


“  *(1> 

{2) 

A 

1.10 

3^54 

0.9* 

OjUJ  j 

6495 

0.91 

>9230 

0.17 

1547 

1.06 

Out 

3:!32 

0.93 

0  89 

1441 

1.21 

0 • 

t 

>1004 

1.20 

Table  1  The  A  Value  Obtained  by  Numerical  Integration 

Key:  1.  Specimen 

2.  (Microstrain) 
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MKO(Bt)  (!) 


Fig.  7 


The  4> ( 9 )  Function  of  Sample  Q 
Different  Load  Levels. 


Key :  (?,)Angle  0  (degrees) 


IV.  Analysis  and  Discussion 

(1)  In  the  direct  ratio  of  the  notch  root's  maximum  deform¬ 
ation  work  density  WQ  to  parameter  J/ p  (formula  8) ,  the  propor¬ 
tionality  coefficient  A  is  approximately  1,  for  a  specimen  of 
a/p==2,A,555  l  and  for  a  specimen  of  a/p  =4. 4,  A  22  1.2. 

Specimen  Q20b  *s  ^irec^ly  clamped  but  because  it  is  not 

clamped  firmly  the  specimen  slips  in  the  clamp  hold.  When 

measured  from  the  auxiliary  strain  piece,  we  can  see  that  the 

specimen  endures  a  relatively  large  bending  moment.  Even  though 

this  be  the  case,  the  (9)  and  A  values  of  the  Q»rt.  are  close  to 

20b 

those  of  the  Q2Qa' 

Besides  this,  tests  initially  proved  that  the  A  value  tends 
to  increase  following  the  increases  of  the  a/p.  Further,  we  can 
infer  that  A  is  related  to  the  hardening  characteristics  of  the 
material.  Yet,  when  applied  in  engineering,  A  can  be  viewed  as 
a  constant  approximately  equal  to  1 . 

(2)  When  the  specimen  notch  is  still  in  a  small  range  yield 
condition,  we  can  use  the  calculation  of  the  stress  strength 
factor  to  obtain  the  J-integral  value.  Thus,  when  establishing 

an  external  load,  the  specimen  and  notch  geometrical  size  are 
related  to  the  notch  root's  maximum  deformation  work  density. 

If  we  assume  that  this  experimental  calibration  is  in  a 
plane  stress  situation,  then  after  using  Weiss's  correction  of 
the  plastic  zone  in  the  calculation  of  we  can  obtain  con- 
sistant  results  (see  table  2) . 
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*  # 

m 

2) 

*0 

(3)  («&$> 

(4)  <»if/a**> 

/• 

J 

l.ll 

89? 

0.0028 

0.054 

0.056 

2.22 

1869 

0.011 

0.22 

0.22 

3.33 

3554 

0.029 

0.60 

0.58 

4.44 

6495 

0.063 

1.31 

1.26 

5.11 

9230 

0.097 

1.97 

1.94 

5.58 

—11500 

0.128 

2.54 

2.56 

1.11 

1443 

0.007 

0.054 

0.076 

2.22 

3004 

0.023 

0.024 

0.25 

3.11 

5119 

0.047 

0.53 

0.50 

0. 

3.56 

6437 

0.062 

0.74 

0.66 

4.00 

8334 

0.037 

1.01 

0.94 

4.22 

9424 

0.100 

1.16 

1.08 

4.67 

—11400 

0.127 

1.54 

1.37 

Table  2. 


Key: 


1 .  Specimen 

2.  (kilograms/millimeter  ) 

3.  (Microstrain)  _ 

4.  (Kilograms/millimeter  ) 

5.  (Kilograms/millimeter) 


J  = 


£  £’  l 


a  + 


•[(-£/- ')i 


(10) 


In  the  formula: 

Y  is  the  geometrical  form  factor  of  the  single  edge  crack 
drawn  specimen  stress  strength  factor; 

is  the  maximum  theoretically  elastic  stress  concentra¬ 
tion  stress  dm=Ktcfg  of  the  notch  root; 

K*.  is  the  actually  measured  elastic  stress  concentration 
coefficient ; 

cr  is  the  distant  mean  stress. 


80 


In  this  paper,  for  C^Qa'  Kt=5*66;  and  for  Qg,  Kt=8.18. 

By  substituting  this  into  formula  (10),  we  can  obtain  the  J-inte- 
gral  value  recorded  as  J*  (table  2) .  At  the  same  time,  from 
formula  (8)  we  can  obtain  the  J-integral  value  recorded  as  J 
(table  2)  when  A=1.0  is  used  for  Q2Qa  an<^  A=l*2  is  used  for  0g* 


Actually,  experimental  calibration  is  not  an  ideal  plane 
stress  condition  as  the  Wq  in  J*=ApWQ  should  be  f dd£  . 

In  this  formula,  f (m)  is  in  the  area  of  1.0-1. 2.  Therefore: 


It  is  also-  only  an  approximate  engineering  calculation  method. 
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J-INTEGRAL  EXPERIMENTAL  CALIBRATION  OF 
SHEET  SPECIMENS  WITH  SINGLE  EDGE  NOTCH 


Luo  Ligeng,  Chen  Xianxi  and  Cai  Qigung 
(Central  Iron  and  Steel  Research  Institute) 
Zheng  Minzhung 

(Aircraft  Strength  Research  Institute) 


Abstract 


J-integral  of  sheet  notched  specimen  has  the  simplest  form  when  the  no¬ 
tched  surface  has  been  taken  as  the  integral  contour,  i.  e. 

/  -  J  Wdy=  W(  0  )  P  co&QdQ 

it  can  be  still  simplified 


J  ^APIV, 

where 

JS/2 

o  <P  (  0  )dQ 

and 

<P(8)~  -cos  8 

it  can  be  seen  that  the  ratio  of  J-integral  to  radius  of  curvature  P  is  propor¬ 
tional  to  the  maximum  deformation  work  density  W0  at  notch  root. 

The  specimens  were  cut  from  hot  rolled  sheet  of  aluminium  alloy  LY12 
with  yield  strength  °o.j=a11.2  kgf/mm1.  The  specimen  sizes  are«  thickness  B  — 
30mm.  width  150mm.  notch  depth  a=40ram,  and  the  radiuses  of  not  h 
curvature  P^^Omm  and  9mm.  respectively.  The  resistance  strain  gauges  of  1 
by  1mm  were  adhibited  on  the  inner  surface  of  notch  at  angles  9  being  O', 
10  ,  20  ,  30’,  45‘,  60",  75"  and  90’,  respectively.  And  then,  the  measure¬ 
ment  of  strain  distribution  s(  9  )  at  various  load  levels  was  made  and  defor¬ 
mation  work  density  distribution  W'(Q)  was  derived  from  stress-strain  curve. 
The  value  A,  at  last,  has  been  determined  by  means  of  numerical  integration 
and  the  approximation  function,  and  varied  in  a  range  of  0.9  to  1.2. 

Owing  to  the  limitation  of  the  experimental  accuracy  and  the  difficulty 
in  determining  exact  degree  of  plane  strain,  the  agreement  of  experimental  ca¬ 
libration  of  J-integral  in  the  present  work  with  J  converted  from  equivalent 
stress  intensity  factor  K,lt  considering  Weiss’s  correction  of  plastic  zone  can 
be  considered  as  an  approximation  for  engineering. 


ON  THE  DESIGN  OF  TRANSONIC  TURBINE  CASCADE  BY  HODOGRAPH  METHOD 
Ling  Zhlguang 

(Institute  of  Engineering  Thermophysics  Academia  Sinica) 

Xin  Shaokang  and  Zhu  Shican 
(Fudan  University) 

Abstract 

The  design  and  realization  of  a  transonic  turbine  cascade 
witn  a  low  energy  loss  coefficient  is  one  of  the  crucial  pro¬ 
blems  for  raising  the  performance  of  transonic  turbines.  This 
paper  deals  with  the  inverse  design  problem  of  the  transonic 
turbine  cascade  in  the  hodograph  plane  by  the  finite  area 
method.  In  the  first  part  of  this  paper,  the  governing  equa¬ 
tion  is  transformed  into  a  symmetrical  form  which  can  be  adopted 
to  obtain  the  upstream  singular  solution.  The  boundary  condi¬ 
tion  connected  with  the  existing  discontinuity  in  the  hodo¬ 
graph  plane  is  also  analyzed  and  afterwards  verified  by  num¬ 
erical  tests.  The  next  part  describes  in  detail  the  method  of 
solving  the  stream  function  field  by  means  of  the  finite  area 
method,  including  the  integral  transformation  of  the  governing 
equation,  the  finite  area  scheme  and  the  particular  corner 
finite  area  scheme.  The  accuracy  of  the  solution  is  also  briefly 
analyzed.  The  worked  out  equation  is  quite  simple  and  quick¬ 
acting.  Calculation  results  prove  that  the  solution  is  quite 
stable  and  it  is  not  necessary  to  modify  the  mesh  near  the  sonic 
line  or  other  places.  Calculations  also  show  that  the  choice  of 
the  position  of  the  critical  point  on  the  sonic  line  has  a 
definite  effect  on  the  flow  field  and  the  blade  type  molded  line 

‘Trans  1  at  or  *  s  Mote:  Subscript  oo  should  be  read  as  ®. 
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form  should  be  given  careful  attention  in  design. 

I.  Preface 

The  transonic  turbine  has  great  potential  for  raising 
efficiency  and  one  of  the  crucial  problems  is  the  molding  and 
designing  of  a  highly  efficient  transonic  turbine  cascade.  Aero¬ 
dynamic  flow  around  body  calculation  methods  such  as  the  time 
correlatio  n method  and  relaxation  method  are  all  forward  solu¬ 
tions  for  various  existing  transonic  turbine  cascades.  In  order 
to  attain  a  transonic  cascade  with  a  low  energy  loss  coefficient, 
if  there  are  only  repeated  improvements  of  the  blade  type  molded 
lines  based  on  aerodynamic  positive  calculations,  the  goal  will 
still  not  necessarily  be  able  to  be  attained.  This  is  because 
the  sonic  line  position  and  its  front  and  rear  flow  fields  are 
all  able  to  directly  affect  the  efficiency  and  performance  of 
the  shock  wave  component  and  transonic  cascade.  Sometimes  the 
application  of  artificial  viscosity  methods  such  as  the  time- 
correlation  method  are  considered  for  the  calculation  of  shock 
wave  discontinuity  and  the  resulting  shock  wave  position  is  a 
region  from  which  it  is  difficult  to  more  accurately  determine 
a  corrected  and  improved  direction.  The  test  research  initially 
revealed:  if  we  could  cause  uniformity  in  front  and  behind 
the  transonic  turbine  cascade  sonic  line,  especially  in  the 
throat  section  flow,  rationally  control  the  speed  gradient  and 
outlet  flow  field  and  rationally  design  a  molding  line,  then 
it  is  possible  to  decrease  the  viscosity  shock  wave  loss, 
realize  a  weak  shock  wave  and  even  a  transonic  wing  and  cascade 
without  shock  wave[l -3  ]. Because  of  this,  we  have  reason  to  be¬ 
lieve  that  the  transonic  turbine  cascade  is  not  very  appropriate 
for  arbitrary  analytical  curve  modelling  resembling  the  subsonic 
turbine  cascade  but  should  adopt  a  general  position  in  accord¬ 
ance  with  the  predetermined  profile  flow  velocity  distribution 
and  sonic  line  and  use  the  inverse  problem  method  for  the  mold¬ 
ing  design.  After  obtaining  the  initial  molded  lines,  the 
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direct  problem  method  is  then  used  for  further  computation 
checking.  This  type  of  inverse  problem  is  carried  out  repeat¬ 
edly  and  hopefully  obtains  a  highly  efficient  transonic  tur¬ 
bine  cascade  series. 

The  hydrograph  method  takes  the  speed  component  or  speed 
molded  length  and  velocity  vector  direction  as  the  independent 
variable  whereby  we  obtain  the  dominant  equation  of  a  linear 
partial  differential  equation  which  simplifies  the  solution. 

There  are  already  many  results[4-7]attained  on  this  aspect  of  the 
transonic  wing.  In  recent  years,  G.Karadimas  and  D.E.  Hobsont8! 
have  carried  out  effective  work  in  extending  the  above  results 
to  the  transonic  cascade. 

As  for  the  inverse  problem,  the  method  provided  for  the  pre¬ 
determined  profile  velocity  is  directly  related  to  the  equation 
and  solution.  If  we  use  the  velocity  distribution  on  the  given 
pressure  surface  and  suction  surface  along  the  x  axis:  71 p= 
f  (x)  and  A  =f  (x) ,  A  is  the  dimensionless  velocity  V/A  _  . 

There  is  also  a  given  distribution  along  the  arc  length  direction. 
This  type  of  supply  method  is  relatively  intuitive  and  seems  to 
be  advantageous  .  to  the  formulation  of  the  control  and  load  of 
the  boundary  layer  loss.  Yet,  this  is  also  conditioned  because 
the  profile  form  and  curvature  are  unknown.  Therefore,  it  seems 
to  be  insufficient  to  guarantee  that  the  gas  emission  angle  for 
the  transonic  turbine  cascade  be  able  to  regulate  the  sonic 
line  location  form  and  a  fixed  level.  The  hodograph  method 
takes  V  and  the  velocity  vector  direction  as  the  independent 
variables  and  supplies  the  relationship  of  the  profile  velocity 
and  profile  inclination.  Aside  from  simplifying  the  equation  and 
the  other  mentioned  advantages,  when  seeking  the  solution,  the 
sonic  line  acts  as  one  of  the  supplied  boundaries.  Thus,  we 
could  also  prepare  beforehand  a  fixed  control  for  the  sonic  line 
location  form.  After  gaining  experience,  we  could  make  some 
restraints  beforehand  to  guarantee  the  gas  emission  angle. 
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Yet,  the  supplied  V  and  6  distribution  on  the  velocity 
surface  form  a  solution  region  and  its  geometrical  shape  has  a 
certain  arbitrariness  and  complexity.  There  are  two  types  of 
methods  used  for  this:  one  is  to  transform  this  solution  region 
to  another  conversion  velocity  surface  causing  it  to  become  a 
rectangular  region  (see  appendix  1) .  The  other  is  even  more 
effective  and  uses  the  finite  volume  method  (binary  is  the 
finite  area  method).  At  present,  the  finite  volume  method  has 
attained  to  very  good  application  for  solving  the  wing  and 
flow  around  the  cascade.  For  example,  Jameson  used  it  to  solve 
the  wing's  transonic  potential  flow[ll]and  McDonald  and  Denton 
separately  used  it  to  solve  the  time-correlation  relation  of  the 
cascaded  12-0}  The  finite  volume  method  showed  the  conservation 
equation  as  an  integration  pattern  causing  the  form  to  be  con¬ 
cise  and  thus  raising  the  stability  of  the  solution.  At  the 
same  time,  a  small  volume  unit  component  mesh  could  be  conven¬ 
iently  deployed  and  divided  in  a  rational  smooth  solution  area 
wherein  we  obtained  a  relatively  good  discrete  approximation 
which  was  advantageous  to  the  calculation  of  the  complex  geo¬ 
metric  shape.  Because  of  this,  we  applied  the  finite  area  method 
to  the  transonic  turbine  cascade  and  even  to  solving  the  in¬ 
verse  problem  of  the  sonic  line's  velocity  surface. 

II.  Dominant  Equation  and  Boundary  Conditions 

1 .  Dominant  Equation 

If  the  shock  wave  in  the  cascade  is  very  weak,  the  effect 
on  the  boundary  layer  is  also  very  small.  We  can  write  mass 
discrete  and  irrotational  equations  for  the  plane,  irrotational 
and  permanent  flows: 


-3J-  (p^-)+-Jr(pF')=  0 


_i_  F.-—  Fy=  0 
ay  w 


Now,  the  independent  variable  changes  to  V  and  0  and  introduces 
stream  function  ^  and  potential  function  $  : 


d<&=V  ,dx+V  ydy 
,  d^= .dy 


Held  on  the  physical  surface: 


d*=  -y-(d<P  +  ) 


After  arrangement,  the  separate  derivations  of  0  and  V  on  the 
velocity  surface  are: 


acj>  _  V  a^_ 

ae  P  »V 


dV  ~  p 


•*»>***“  ■■■ 


After  merging  and  becoming  nondimens ionalized,  we  can  obtain; 


f i 


In  the  formula 


<?- 


J.*  _  n _ 

v  Pm^Jcos0m  '  r  i 

A  o  _  *  ~  1  . 

(  1  -BA*)*  '  "  k  +  i »  c 


_J  -A* 

~( i  -  /?/.*)"•'  * 

i  (1) 

=  t  %mmo 


Key:  1.  "fc  is  the  blade  pitch 

Formula  (5)  is  the  linear  partial  differential  equation  needed 
for  the  solution. 

Formula  (5)  can  be  further  transformed  into  a  symmetrical 
form.  Letting  £  =P,,r$  and  jj’p _  (1-B  A  =  f  (  A.  )  ,  after 

using  the  new  equivalent  independent  variable  7\J  to  substitute  for 
A.  ,  we  obtain: 


^  d// 

(  c  ) 

* 

i 

k 

■Ife 

( 7 ; 
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In  the  formula 
which  is 


and 


X' 


■f 


i/l  —  M 


BX, 


1+B,  17 

log[| 


y/  1  —  )  s  ~ "j/^1  —  BX* 

Vi  -x2+K/i-5xjf 


>/t 


fiy/1  —  X*— •»/  /?  C 1  -  BX*) 
fl/l 


_j 
2  ✓  C 


+  C 


Thus,  formula  (5)  becomes: 


1(4*) 


(  8  ) 


2.  Boundary  Conditions 

After  utilizing  the  stream  function  from  figure  1,  the 
boundary  conditions  of  y *  =0  on  the  suction  surface  along  the 

blade  surface,  V  *=1  on  the  pressure  surface  and  i|>  =f(0)  on  the 

P  c 

sonic  line  can  be  determined  by  using  the  singular  solution  of 

a  mixed  type  equation [6]. 


Fig.  1. 


Fig.  1. 


Key:  1.  Velocity  surface 

2.  Physical  surface 

3.  Sonic  line 


Up  to  the  present,  the  difference  between  the  various  singular 
solutions  lie  in  the  accuracy  of  their  compression  functions. 
Other  approximate  hypotheses  were  all  the  same  but  in  the  area 
of  the  sonic  line  the  differences  of  the  compressibility  func¬ 
tions  were  very  small.  When  we  used  Germain' s[10|  solution,  which 
shows  that  there  is  a  critical  single  point  6c  nozzle  on  the 
throat  part,  then 


'  (i3,-eJls-(8,-9.),/s 


The  upstream  boundary  conditions  are  the  distant  overtaking 

flows  A-qo  and  Qqo  which  forms  single  point  I  on  the  velocity 

surface  from  the  cycling  of  the  flow.  Further,  after  selecting 

the  blade's  "geometric"  stagnation  points  A.^  and  A  ,  we  can  then 

determine  the  solution  area  IM.  A,  C.  C  AMI  on  the  velocity 

liioooo 

surface.  To  sum  up,  the  problem  is  the  solution  of  "degenerated" 
elliptic  equation  (5)F*N*^  or  equation  (8). 

Because  the  equation  is  linear,  its  solution  can  be  obtained 
from  a  combination  of  the  regular  solution  and  singular  solution, 
which  is  y  =  +  "V  .  In  this  equation,  "’p/  is  the  regular 

solution  and  ^  is  the  singular  solution  near  the  upstream  I 
point.  Therefore,  the  solution  equation  changes  to: 

F.N.l:  Because  it  is  necessary  to  directly  solve  the  singularity 
sonic  line. 
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L  (>!>' )  =  /.  (H>*)  -  L  ( H'- )  =  -  L  (Hr.) 


(9) 


we  used  an  approximation  of 

A 


7~ 

oo 


W= 

^  oo 

formula  (8)  changed  in  the  Laplace  equation 
directly  attain  the  Fenain  results: 


in  formula  (8)  and 
From  this, we  could 


4-= 


(0 

2* 


tgfl- 

‘  2**/  i  -  Mi 


log/?  +  C 


(10) 


In  the  formula 


/?  = 


+  l-Xi"’ 


/ 


n  -n„  ^ 

✓  l  ->.i  J 


Now,  the  solution  area  is  A.C.C  A  M  M,A,  and  the  boundary 

^ _ ^  lloooll 

conditions  are:  7= for  A_C  ;  'V  =1-  ^ 

OO  a  o  *  t 


o  o' 


for  A, C, ;  and 
oo  1  1 


i'*  t  -  y  for  £c\  is  necessary  to  further  discuss  the 
c  oo  1  o 

A-M.M  A  line  because  the  M.M  area  function  is  non-continuous. 

1  l  o  o  1  o 

It  jumps  from  the  'f  =  1  in  the  area  to  the  'V  =0  in  the  Mq 
area  and  is  also  on  the  A  =0  line:  for  A.M.  r  *  1  -tfi 

OO  11  r  ' oo 


^or  AoMo  In  order  to  maintain  a  continuous  ^  on  the  solution 

area  boundary  it  is  necessary  to  establish  jIM,  =tIm  .  This 

i  o 
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can  be  realized.  For  example,  using  M.  =  =0  it  is  then 

I  O  / 

necessary  that:  ^  lM-=l  and  To  do  this,  it  is 

OO  1  oo  o 

only  necessary  to  use  the  C  value  in  formula  (10) : 

C=  1  — y )  co<co4 

<2) 

C - £  ®>co^ 


Key:  1.  When 

2 .  When 

Then  the  goal  can  be  reached. 


Fig.  2. 

We  can  also  use  reference  [1]  which  takes  Wm  =  "fl M  =1.  At 
,  i  o  l 

this  time:  f  M, =0  and  y  I  M  =1  and  it  is  necessary  to  use 
oo  l  '  oo  o 

c=-' 27f  when  w  a  and  C=”  2tT  -1  when  w  ^“a* 

Aside  from  this,  we  can  also  use  T/M  =  y  /  M.  _ _  _ 

o  1  as  another 


arbitrary  value  such  as  0.5  etc. (see  fig.  2).  From  this  comes 

the  corresponding  C  value  and  given  ^A-  Numerical  tests  proved 

that  no  matter  what  the  value  of  YIm  =  "'W  M.  ,  the  results  were 

o  1 

the  same.  We  could  also  understand  why  the  final  result  was 

y*.  1  .  The  above  result  was  obtained  assuming  IM  was 

oo 

obtained  in  a  straight  line.  Actually,  IM  represents  the  upstream 
flow  field  and  any  &  origin  can  be  selected  along  IM. 

Ill .  Solution  Method 

We  used  the  finite  area  method  to  solve  the  boundary  value 
problem  of  formula  (9) . 

1 .  Transformation  of  Equation 

In  order  to  be  appropriate  for  solving  the  finite  area 
method,  the  integral  of  formula  (9)  was: 


III  k  (Plf,i)  +  ]  <f0dX=JI/l(  9  ’  x  )</0dX 

D  D 


In  the  formula,  f  ,  (-0j  ,  A)  =L ( 'f  ).  After  applying  the  Green  form- 

1  oo 

ula : 

(P^d\-QYidQ)  =<fi  +  l)d$d\  (11) 

J  dD  ~  Q 

In  order  to  simplify  calculations  below,  Y'  is  written  as  Y . 
2.  Treatment  of  Format 

Using  a  non-equidistant  A  =const  line  (actually  it  is  not 
absolutely  necessary  to  use  an  equal  A  line,  yet  because  A=0 
and  A  =1  are  boundaries,  the  use  of  an  equal  A  line  is  more  con¬ 
venient)  and  an  equally  separated  point  joining  line  along  the 
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A^const  line,  we  cut  up  the  solution  area  into  an  arbitrary 
tetrahedron.  The  stream  function  t  was  defined  in  each  format 
center  area.  For  each  tetrahedron  mesh  in  the  solution  area, 
peripheral  integrals  in  formula  (11)  can  be  changed  into  line 
integral  sums  on  the  four  sides. 


i 

the 


Fig.  3. 

For  example,  the  2-3  side  in  figure  3  has: 

Pif»  j:  n„ds  +  Qy,  J  ^  n  ds 
”  X*)  +  — 9j) 


In  the  formula,  P,  Q,  f  Q  and  f  are  all  mean  values  in  the 
integration  area.  We  took  point  P  and  Q  values  in  the  area  as 
P  and  Q.  Because  "f  was  only  defined  in  the  format  center  areas, 
it  was  necessary  to  provide  the  expression  of  derivatives  f  0 
"f  y  Using  the  directional  derivative: 


at) 

dl  ~  *0  di  il 


tMI+wk 


If  for  point  ^  of  the  2-3  side  in  fig.  3  we  take  the  1  direction  as 
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the  ee1  direction  and  2-3  direction,  then: 


*  ^,J'C0SU‘ +  +«*"«. 


-  v7^J“”^,,;C08P'  +  'I’ktSinP, 

£\«3j 


In  the  formula, 
points; and  &r 


S^2  is  the  distance  between  the  two  adjoining 


ee 


1 

format  centers.  VJvQj  is  the  point  area  7a,  and  *<S> 
point  (^  area  %  (see  fig.  4).  when  a^=0  (using  the  equal  A  line 
to  divide  the  mesh) ,  we  can  obtain 


is  the  distance  between  the  two  adjoining 


Ta.  and  "'La  is  the 


'll*  D 5 


l|ln  — ' 


^00  = 


_ 


Arsing, 


ctgPi'Ju® 


(12) 


Fig.  4. 

Taking  +,-  *  (  t.*  t.,+  *..♦  %  >  and 


4A»„sinP, 


r»l 


:tgPi 


Therefore,  the  2-3  side  can  be  written  as: 


/«-  P  (!.)-%  ■**-(Xt-X,)  +  0  (X,) 


— CtgPj 


Ar##j 


](9i-9s) 


4A«uSin0i 


Fere,  ^3'  93 '  ^2  an<^  e2 are  the  A 
and  P  (  7l-e)  and  Q  (  A.  are  the 
area . 


and  ©  values  of  points  3  and  2, 
P  and  Q  values  of  the  7L  =  A 

e 


Similarly,  point  (?)  of  the  1-2  side  has: 


A*.. 


4As,, 


— +T*«« _ ctgp4 

T‘*  Af  oinfl  cl6l-'4  .  A. 


Ar„4sin04 


4Aslt 


Because  7^  =  ^2'  therefore 
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b  and  x®  are  analagous.  When  the  four  sides  are  added  to 
gether  we  obtain: 


+'», — / 1  (/•»*  0.) 


_Afi.  t_A«u 
2 


ijr.smj}, 


In  tho  formula,  a,  is  the  coefficient  of  V7  and  r 

1 ;  e  e ; 


The  boundary  infinitesimal  has  eight  types  of  different 
situations  (Fig.  5)  wherein  it  is  necessary  to  treat  and  write 
out  their  integration  formulas.  For  example,  for  situation  (p, 
from  Fig.  6,  ‘'J' =1-  on  the  1-4  side,  and  V  =  fc  (©)-  "*P00 

on  the  3-4  side.  Therefore 


Fig.  6. 


/,-Q(l)(0,-04) 


^l>.(0j) o  4«(0») ~4i-t~ ( 1  - 

"  4*“ 


and 


'l'«x= 


'I'.  -  (  1  -  '!>-*) 

^  A^si  4-  j 


l(’-  r  —  'i’-X 
As4lsi«03 


“CtgPj 


i|>,—  (j  —  ^-t) 

'(*!*) 


Therefore 


/»=  P  (X.) 


'!).-(  l  -J'-i) 
^  Ags«  +  AJn^ 


(X,-X4)  +  9  (X.) 


A*41sint) 


Finally,  the  related  simultaneity  of  the  infinitesimal  mesh 
defines  the  linear  algebraic  system  of  first  order  equations 
in  each  format  center  1*  and  then  the  solution  is  carried  out. 

3.  Annotation  of  Several  Points 

(1)  For  the  most  part,  the  finite  area  method  is  carried 
out  for  a  system  of  first  order  equations.  For  example:  equa¬ 
tion  (5)  can  also  be  written  to  form  a  system  of  first  order 
equations 


1  d$>  _ 

Q  *9 

30  ■  n 


=■  0 


s  0 


After  integration 


(13) 


At  this  time,  the  mean  values  of  the  two  adjoining  mesh  format 
centers  are  used  for  the  value  on  the  infinitesimal  area  bound¬ 
ary.  The  presently  used  equation  (5)  has  direct  integration  and 
utilizes  the  derivative  value  on  the  boundary.  The  coordinates 
of  the  four  angle  nodal  points  use  the  mean  values  of  the  four 
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mesh  format  center  .  The  accuracy  is  raised  and  some  are  similar 
to  the  Lax  format.  The  use  of  the  4  point  means  is  also  advant¬ 
ageous  to  stability.  Aside  from  this,  the  use  of  a  system  of 
first  order  equations  can  also  save  on  the  machine's  storage 
quantity. 

(2)  In  the  solution,  the  first  order  derivative  approx¬ 

imation  uses  the  difference  quotient  alternate  approximation  of 
the  1  and  Q-directions  and  the  accuracy  is  in  the  first  to 
second  order  range.  When  the  mesh  degenerates  to  a  square  it 
has  second  order  accuracy  equal  to  the  central  difference. 

From  Fig.  3,  when  B^=90°, 


an  !  i ,  j  *  1/2  AX. 

Hi  1  ^4>ni./-ik# 

,■  *i/i,  j  A0 


(•3)  We  can  also  use  the  weighted  residue  method  of  taking 
the  weight  function  as  1  for  the  integration  of  formula  (8) . 

If  the  weight  function  acts  as  the  form  function,  mathematic¬ 
ally  speaking,  this  corresponds  to  the  finite  element  method. 

(4)  is  the  singular  solution  "near"  the  upstream 

infinity  area. .We  continued  the  use  of  the  method  in  refer¬ 
ence  [1)  and  employed  it  for  the  whole  solution  area;  although 
this  method  was  much  more  convenient  to  use  than  the  hyper¬ 
geometric  function  for  seeking  the  effects  of  the  upstream 
single  point,  yet  the  approximation  was  larger.  Therefore,  in 
the  next  step,  we  calculated  the  problem  of  the  first  category 
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boundary  value  conditions  in  carrying  out  a  direct  solution  of 
formula  (8)  or  (5)  on  a  fractured  surface. 

IV.  Calculations  and  Discussion 

Based  on  the  above  mentioned  method  worked  out  on  a  computer 
program,  we  obtained  the  ^  field  in  the  solution's  linear 
algebraic  series  of  equations.  After  obtain  the  Y  field  we 
used  numerical  integration  to  directly  calculate  the  s  surface, 
p  surface,  and  x  and  y  sonic  line  coordinates.  The  coordinate 
origin  was  taken  on  the  intersecting  point  of  the  s  surface  and 
sonic  line.  For  the  supersonic  part  we  continued  the  use  of  the 
characteristic  line  solution  method. 


Fig.  7(a)  The  Relationship  of  the  A -8  Distribution  System 


Ill 


/e<  =  6a.8 
\  *  70.0 


Fig.  7(b)  Profile  Calculation  Example 

Fig.  7  is  the  calculation  results.  The  cascade  nodal  dis¬ 
tance  t is  taken  as  1,  the  upstream  overtaking  flow  is  M  =0.15, 

o 

the  axis  gas  intake  is  aQ=90°  and  the  gas  emission  angle  is 
a^  q  15°.  Fig.  7(a)  gives  the  _9  distribution  and  Fig.  7(b) 
is  the  sought  transonic  turbine  guide  blade  profile. 


(1)  In  the  relationship  of  the  given  velocity  distributions, 

at  the  sonic  line  A.  =1  area,  G  =72.3°  and  9  =68.7°.  The  calcula- 

P  s 

tions  practically  proved  that  the  selection  of  the  critical 


single  point  ©c  was  not  only  related  to  the  performance  of 
the  supersonic  part  form  but  the  effect  on  the  profile  of  the 
front  surface  of  the  sonic  line  and  the  form  and  position  of 
the  sonic  line  was  rather  great  and  quite  sensitive.  Fig.  7(b) 
shows  the  results  created  by  three  types  of  0c  values  wherein 
0c=68.8°  is  taken  as  relatively  good.  It  can  be  seen  from  this 
that  when  we  want  to  obtain  a  uniform  flow  in  the  throat  area, 
the  subsonic  line  form  is  very  important.  After  providing  pres¬ 
sure  surface  p  and  the  A relationship  of  the  suction  surface, 
following  adjustment  the  ©c  can  obtain  a  relatively  good  value. 
Concrete  calculation  examples  have  also  showed  that  in  a  certain 
overtaking  flow  M  number  sphere,  the  9C  point  near  the  back  arc 
will  cause  the  sonic  line  to  become  even  more  level  and  straight. 
This  is  equal  to  causing  the  exit  back  arc  to  approach  a  straight 
line  shape. 

(2)  After  the  ^  *  field  is  known,  the  return  physical  sur¬ 
faces  x  and  y  require  separate  integration  formulas  for  the 
p  and  s  surfaces: 

{[P  +  Q{-w)  005  8  (14> 


(15) 


For  the  sonic  line 


[-^-cosO  -^-sin0  ]rf8  (16) 

J  ^  [-^-ain  0  +-^-cos  0  j<f0  (17) 

_  Poo 

In  the  formulas,  the  C=  P  A  ooC0S®oo  is  obtained  from  the 

given  overtaking  flow  conditions.  In  order  to  satisfy  the  profile 
seal  condition  there  should  be 

(*■ - *,)  | ,  +  (xm - xe )  |, =  (xm - X,  „)  | , 

In  the  formula, n  indicates  the  integration  end  points  of  the 
p  and  s  surfaces,  m  is  the  end  point  of  the  sonic  line,  x 

op 

and  yQp  are  the  coordinates  of  the  p  surface  integration  begin¬ 
ning  points  and  x  =x  .  If  i  is  used  to  indicate  the  flow 

op  mo 

points  on  the  p  or  s  surfaces,  then  the  blade  thickness  distri¬ 
bution  can  be  shown  as : 

Si-  t  +iyi-y^\,-(yi-y„)\,-(.ym-yt)\. 


During  the  actual  calculation  procedure,  it  was  necessary 
to  repeatedly  adjust  and  modify  the  original  data  (the  7V  -8 
distribution  of  pressure  surface  p  and  suction  surface  s)  or 
consider  the  relationship  of  the  profile's  form  velocity  and 
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inclination  changes  so  as  to  obtain  a  relatively  good  seal  and 
to  satisfy  the  thickness  distribution  of  the  profile.  After 
gaining  a  certain  amount  of  experience,  this  was  not  difficult. 
Because  the  function  P  value  had  singularity,  when  there  was 
actual  integration,  we  could  try  to  shift  to  the  upper  half  of 
the  mesh.  Although  this  caused  the  leading  edge  to  have  a  small 
number  of  breaks,  this  problem  did  not  have  a  very  large  effect 
when  there  were  many  mesh  nodal  points. 

(3)  When  using  the  finite  area  method  for  solution,  it 
was  necessary  to  divide  the  mesh  for  the  solution  area  and 
solve  the  mesh's  central  values.  The  mesh  was  divided  in  detail 
and  solution  accuracy  was  raised.  As  regards  the  solution  re¬ 
sults,  if  a  relatively  rough  mesh  was  used,  the  results  were 
completely  stable  near  the  A.  =1  area.  The  same  holds  true 
for  the  A.  =0  area.  In  this  way,  it  was  not  necessary  to  have 
relaxat iv difference  solution  in  some  places.  For  example, 
with  the  us<5  of  under-meshing  near  the  sonic  line,  the  mesh 
was  finer  in  the  A  direction  or  a  mesh  with  an  undetermined 
value  was  used.  This  also  used  the  special  characteristics  of 
the  finite  area  method  and  was  one  of  the  advantages. 

V  Brief  Summary 

(1)  This  paper  presented  the  use  of  the  finite  area  method 
on  the  velocity  surface  for  the  solution  of  the  transonic  tur¬ 
bine  cascade  as  well  as  the  inverse  problem  results  of  the 
sonic  line,  a  concrete  description  of  the  treatment  of  the  equa¬ 
tion,  the  solution  method  steps  and  the  handling  of  the  format. 
It  also  discussed  the  accuracy  of  the  solution.  After  using 
this  method  for  numerical  calculation  and  analysis,  the  program 
was  concise  and  calculation  time  was  short  which  could  be 
appropriate  for  a  complex  area  form.  Moreover,  it  was  not  nec¬ 
essary  to  have  more  meticulous  treatment  near  the  sonic  line. 
This  point  is  an  improvement  over  the  method  found  in 
reference  [ 1 ] . 


L5 


I 


(2)  Calculations  showed  that  the  selection  of  critical 
single  points  of  the  throat's  sonic  line  which  have  a  rather 
great  effect  on  the  sonic  line  form  and  profile  form  should  be 
given  careful  attention  in  modelling  design  calculations. 

Appendix  The  Transformation  of  the  Solution  Area  on  the 
Velocity  Surface 

9 _ 9 

In  taking  i  *  and  n=Y  as  the  transformation  of 

the  coordinates  on  the  velocity  surface  (Fig.  8) ,  we  had: 


a  1  a  a  _  a  ax  a 

a0  =  0,-0,  ax'  a\  ax  aX  ~  ay 


p  X 

Letting  g=  -5*7,  we  obtained  from  the  relationship  in  the  coordin¬ 
ates: 


n  _  (e  -e,)e;-<e  -e.)e; 

9  (0,-0,)*  "■ 


in  the  formula  e'p  and  0's  are  the  0p  and  «s  partial  deriva 
tives  for  A-  .  in  the  same  way,  we  can  obtain: 
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g,  m ^4-  „  Mm, + M--  M- »?) + 8'( " ee' 4  M' + 0,9 "  %) 

+2o;o;(o,-20+8>)+2e;i(e  -&,)+2or-(0 


Letting  gy=  ^  and  Q^=  ^  ,  on  the  converted  velocity  surface, 
equation  (5)  becomes: 
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THREE  DIMENSIONAL  STRESS  ANALYSIS  FOR  A  SHROUDED  AIR-COOLED 
TURBINE  BLADE 


Shung  Changbing  and  Ziao  Junxiang 

(Beijing  Institute  of  Aeronautics  and  Astronautics) 

Abstract 

Up  to  now,  there  are  still  no  accurate  and  reliable  methods 
for  computing  the  stress  field  of  a  shrouded  air-cooled  turbine 
blade  because  the  cooled  blade  with  a  complex  figuration  works 
under  high  temperatures  and  centrifugal  loads;  moreover,  the 
distribution  of  aerodynamic  loads  is  rather  non-uniform.  Gener¬ 
ally,  such  problems  which  occur  in  practice  are  handled  by 
means  of  experiments  of  statistics. 

In  attempting  to  solve  these  problems,  we  applied  the  high 
order  three  dimensional  isoparametric  elements  for  describing  the 
complex  configurations  of  blades.  In  order  to  solve  the  high 
order  linear  system  of  equations,  we  have  used  the  frontal  method 
with  high  accuracy  and  considerable  economy  of  main  memory  space. 
The  high  order  bicubic  Coons  interpolation  function  is  also  used 
to  fit  the  space  curved  surface.  In  order  to  verify  the  present 
method  and  the  reliability  of  the  given  program,  we  have  compared 
the  results  of  numerical  computations  with  those  obtained  from 
experiments  on  the  revolution  of  the  real  blades  at  high  speeds, 
and  the  comparison  confirms  the  adaptability  and  the  satisfactory 
accuracy  of  the  present  method. 


This  method  is  appropriate  for  various  blades  such  as  solid 


or  hollow  blades,  shrouded  or  unshrouded  blades,  blades  with  or 
without  a  shank  root,  and  at  high  or  normal  temperatures. 

I.  Preface 

The  analysis  and  calculation  of  the  high  temperature  fields 
and  stress  fields  of  shrouded  air-cooled  hollow  blades  is  a 
major  research  topic  brought  forth  by  related  departments  in 
China.  When  used  by  an  engine,  it  can  improve  the  performance 
parameter  greatly  and  raise  its  lifespan.  At  present,  there  are 
already  many  types  of  advanced  engines  making  use  of  it.  After 
the  appearance  of  the  finite  element  method,  references  [1]  and 
[4]  had  calculated  the  temperature  field  and  thermal  stress  for 
the  shrouded  hollow  blade.  Yet,  the  three  dimensional  blade  was 
simplified  into  the  handling  of  the  problem  of  the  plane.  Natur¬ 
ally,  there  is  a  great  difference  between  this  and  the  real 
situation . 

A  few  of  the  major  characteristics  of  this  type  of  blade 
are:  firstly,  it  is  composed  of  an  irregularly  shaped  shroud, 
possesses  a  special  profile  and  a  warped  blade  body  as  well  as 
various  special  types  of  cooling  holes  and  shank  roots.  Because 
of  this,  is  is  very  difficult  to  disperse  on  the  curve  side 
element,  to  fit  on  the  special  form  and  to  automatically  pro¬ 
duce  on  the  nodal  point  coordinates.  Secondly,  the  blade  works 
under  a  complex  load.  Its  tangent  velocity  generally  reaches  to 
250-300  meters/second  and  therefore  the  centrifugal  load  is 
very  great.  Moreover,  the  barycentric  shifts  of  each  tangent 
surface  produce  added  moment.  When  the  unevenly  distributed 
aerodynamic  load  is  along  the  arc  direction  and  radial  direc¬ 
tion,  it  is  extremely  difficult  to  simulate  the  boundary  condi¬ 
tions.  The  working  temperature  of  the  air-cooled  blade  (in  front 
of  the  turbine)  reaches  to  over  1500K,  there  is  very  great 
thermal  stress  produced  on  the  side  of  the  cooling  hole  and  the 
high  temperature  area  often  has  burning.  Thirdly,  testing  and 


measuring  are  difficult.  When  the  blade  works  under  high  tempera 
ture  and  high  rotational  speed,  it  is  impossible  to  directly 
measure  the  required  temperature  and  stress  distribution  of  each 
point.  Moreover,  it  is  even  more  difficult  to  simulate  the 
engine's  various  working  conditions  in  a  laboratory.  Fourth, 
theoretical  calculations  are  difficult.  There  are  many  load 
types  that  require  calculation  for  the  blade  and  this  is  very 
difficult.  For  example,  at  present,  the  temperature  boundaries 
and  aerodynamic  boundary  conditions  are  being  studied,  yet  the 
calculation  of  the  temperature  field  and  stress  field  depend  on 
the  accuracy  of  giving  these  conditions.  Furthermore,  the  order 
of  the  series  of  equations  is  also  an  outstanding  problem. 


i’ig.  1.  External  Form  of  the  Shrouded  Air-Cooled  Blade 

Key:  1.  Shroud 

2.  Gas  emission  hole 

3 .  Tangent  plane 
4-  Gas  intake  hole 
5  -  Root 


■  2 

(a)  *#*#<*•  y>  2)»  (6) 


Fig.  2 


Key:  (a)  Integral  coordinate  (x,y,z); 
(b)  Local  coordinate  (|  ,  h  , £  ) 


The  form  function  of  the  20  nodal  points  on  the  indicated 


element  is: 


r 


tfi-jd  +4.)d  +n,)(  1  +to(l.+t«.+t.-  2) 

(  i  =  1,  3,  5,  7,13,  15,  17,  19) 

iV,=  l(  i  - 4‘) ( i  +n0)(i  +£„) 

4 


(  <  -  2,6,  18,  14) 
iV,  =  l(  i  -tl*)(l  +t,)(l  +4,) 

(•=•4,8,  20,  16) 

W#~j(i  -t*)(i  +4,)(i  +n.) 

( «  =  9 , 10,  11,  12) 


(i) 


In  the  formulas ,  %  -  L .  1  ,  n  =  0  .  h  ,  £  =  ^  i  ^  ,  and 

o  l  o  l 

t,  .  ,  Hi  and  are  local  coordinate  values  of  the  i  nodal 

points  - 

We  know  from  the  elastic  theory  that  for  the  geometric 
equation  of  the  three  dimensional  problem  the  matrix  is  used 
to  indicate  it  as : 

{ «  >  =C5) { 5 }  ( 2 ) 


In  the  formula 

Key:  is  the  matrix  of  the  relationship 

of  the  strain  and  displacement; 

ft)  is  the  displacement  vector? 

{£}  is  the  strain  vector 


Its  physical  equation,  when  there  is  initial  strain  [^0] * 


<o>  =  CDH<e>-  U..» 


In  the  formula 


Key:  {<*}  Is  the  stress  vector; 

f£#?  =  aT  {111000>  ; 

T  is  the  temperature; 

a  is  the  material's  linear  expansion 
coefficient ; 

CD)  is  the  elastic  matrix . 


The  element  has  volume  force  and  surface  force  in  simultan¬ 
eous  action  and  the  initial  load  functional  equation  caused  by 
the  temperature  difference  is: 


rr(u,  v,w)=  ?>}*-({ 6 


For  the  whole  area,  then  we  have 


r 

n  - 2.n*=  Y.  4-(<MVc*m>-(<i>r<FrJ 


In  the  formula 


{F),=  {F,y+{F.V+{FAy 

V* 

W-Bftwpxr 

V 

{F.ty~$iNy{<ndv 

A" 


is  the  element's  external  load  (5a) 
vector ; 

is  the  element's  initial  load  (5b) 
vector; 

is  the  element's  volume  force;  (5c) 
is  the  element's  surface  force;  (5d) 
is  the  element's  stiffness  matrix;  (5e) 


... 


rt>  =(rai>  • 


Csl  is  the  nodal  point  shift  vector? 

.  INj  etc. 


When  functional  n  uses  an  extreme  value  of  5u  =0,  then  we  obtain 
the  series  of  equations  required  for  solution. 

III.  Temperature  Field' 

In  order  to  calculate  the  element ' s  thermal  load  and 
thermal  stress,  it  is  necessary  to  calculate  the  three  dimensional 
temperature  field  on  the  blade.  Based  on  the  experimental  data 
found  in  reference  [12]  concerninq  geometrical  correction  along 
the  arc  length  and  mathematical  interpolation  along  the  blade 
height  for  a  certain  type  of  blade  we  obtained  its  boundary 
value.  When  the  boundary  conditions  were  of  the  first  category 
and  third  category,  the  temperature  field's  elliptic  equation  and 
boundary  conditions  were: 

j  +  a'-T  +  a-T  =  Q 

|  ax 1  ay1  as-  (in  volume  V) 


T  z=T„(  x ,  y,  z) 


(on  surface  A^) 


aT 

an 


+  %(x,  y,  z)  =  q  (x,  y,  z) 


(on  surface  A2) 


(  6  ) 


In  the  formula 


A=A1+A2 

T 

j  T  O 
*  n 

nQ  and  q 


is  the  entire  boundary  condition? 
is  the  known  temperature  distribution? 

is  the  external  normal  derivative? 
are  boundary  parameters . 


Using  the  Gauss  divergence  theorem,  we  can  transform  the 
elliptic  equation  into  a  functional  extreme  value  condition: 
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(7) 


»n  =  2 c  (  {5t}  ')tc  k}+ a'd  {bA  r  { yn 

e 

*{&r}r(c/cD<n-{F})=  o 


In  the  formula 

CA'(0*=JjTcC3T(C/ri)TC/r‘CCDc/F  is  the  volume  element  con- 

y.  tribution  to  "total  stiff¬ 

ness"  ; 

CA,1'=Jj  {is/)  {M)r\dA  is  the  surface  element  con- 

tribution  to  "total  stiff¬ 
ness"; 

{F),=  ((iN)qdA  is  the  surface  element  con- 

tribution  to  the  right  end. 


Due  to  the  arbitrariness  of 


formula  (7)  is  equivalent  to 


CfO<T}-<F}  =  0 


(8) 


This  then  is  an  algebraic  system  of  equations  for  the  element's 
nodal  point  temperature. 

IV.  Frontal  Solution  Method 

This  is  a  variation  of  the  Gauss  method  and  its  greatest  ad¬ 
vantages  are  that  its  solution  accuracy  is  high  and  its  time  is 
fixed.  Moreover,  its  programming  is  ingenious  so  that  it  can  de¬ 
crease  internal  memory  and  make  full  use  of  external  memory. 

Its  special  characteristics  in  operation  are  that  it  both  accum¬ 
ulates  and  eliminates,  it  does  not  require  the  formation  of  total 
stiffness  and  the  results  after  elimination  are  sent  to  the  ex¬ 
ternal  memory  and  does  not  occupy  internal  memory  space. 

By  using  fig.  3,  it  is  even  easier  to  explain  its  solution 
process.  In  the  fig.,  when  elements  (JM?  are  already  accumulated. 
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element  ©  is  in  the  process  of  being  carried  out  and  elements 
(§)  -  have  not  yet  accumulated;  nodal  points  with  a  symbol 
have  already  been  eliminated,  nodal  points  with  a  "0"  symbol  are 
in  internal  memory  which  is  equivalent  to  the  tape  width,  and 
nodal  points  with  a  "A "  symbol  have  still  not  begun  action. 
Therefore,  the  activity  variable  "0"  in  the  solution  procedure 
forms  a  "wave  surface"  that  continuously  advances  forward.  This 
is  the  origin  of  the  frontal  method. 


a-  (.'■)  t-.iKfVfflM  r>£  □-  «*£««** 


Fig.  3. 


Schematic  Diagram  of  the  Frontal  Method 

o  is  a  nodal  point  in  the  internal  memory 

x  is  a  nodal  point  already  eliminated  element 

4  is  a  nodal  point  which  still  has  not  begun 

action 

is  an  element  being  produced 
Bl  is  an  element  which  has  already  been  produced 

□  is  an  element  which  has  not  yet  been  produced 


It  is  easy  to  prove  that  the  internal  memory  of  the  frontal 
method  mainly  depends  on  the  width  of  the  wave  surface  and  be¬ 
cause  of  this  it  can  economize  the  internal  memory  better  than 
the  band  width  method  and  divided  piece  tape  width  method. 

Tests  have  proven  its  obvious  advantages  for  high  order  systems 
of  equations  and  high  order  elements.  Its  shortcomings  are  that 
the  program  is  too  complex  and  when  operating,  the  repetition 
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of  internal  and  external  memory  exchange  wastes  time. 


V.  Fitting  and  Automatic  Mesh  Division  of  the  Profile's  Curved 
Surface 

1.  Fitting  of  the  Profile's  Curved  Surface 

The  body  of  the  blade  is  surrounded  by  two  skewed  space  curved 
surfaces.  The  contour  of  the  curved  surface  must  be  accurate 
otherwise  this  will  affect  aerodynamic  performance  and  engine 
efficiency. 

The  divided  three  dimensional  curved  surface  element  must  not 
only  have  the  curved  surface  strictly  coincide  with  the  profile 
but  the  element ' s  nodal  point  coordinate  must  also  accurately 
fall  on  the  curved  surface. 

We  have  test  tried  a  new  method  proposed  in  recent  years  by 
Coons  and  Forest  [13]  using  a  type  of  structured  curved  surface 
for  the  fitting  of  the  space  curved  surface.  If  the  profile  is 
composed  of  many  small  curved  surface  pieces  and  the  interpola¬ 
tion  function's  second  order  derivative  is  continuous  in  each 
of  the  curved  surface  pieces  (as  shown  in  fig.  4) ,  the  domain 
of  the  bivariable  interpolation  function  x(u,w)  is:  0  4  u(  1 

and  0  w  ^  1 .  If  we  now  fix  a  certain  parameter  in  the  x(u,w) 
then  we  obtain  the  common  single  variant  interpolation  function. 


U/+1)  (l+D,  0+i) 


i,i  (»  +  l),/ 


Fig.  4, 


Fig.  4.  Curved  Surface  Pieces 

When  formulating  the  partial  derivative  of  the  second  order 
continuity  in  the  boundary,  we  write  the  matrix  as: 


x(u. 

«')-CF,(a)  G, 

(“) 

F,(  a  ) 

Gt 

(u» 

*(0 

,  0) 

*»(0, 

0) *(0, 

1) 

xm(  0 , 

D 

F0(a»)  \ 

*.(  0 

,  0) 

*«.(  0 , 

0  ):*.(  0  , 

1) 

*«.(  0, 

i) 

<?.(»)  j 

*(  1 

,  0) 

*.(  1 , 

0):*  (  1  , 

1) 

*.(  1 , 

i) 

F,(u») 

*.(  1 

,  0) 

xmm(  1 , 

0  )*.(  1 , 

1) 

*.»(  1  > 

i) 

Gt(w) 

In  the  formula,  the  variable  in  the  0-1  area  is  second  order  con¬ 
tinuity  and  the  related  symbols  are: 

F.Co^F.Cl)-  1 

F0(l)=F1(0)-F;(0)-Fi(l)=F;(0)  =  F,(0)  -  0 

G5(0)=G((1)=  1 

•  Ga(  0  )=(?,,(  1  )  —  C, (  0  )==<?,(  1  )=GI(  1  )=(?!(  0  )  =  0 

x.  =  dx/du  4$ 


Key:  1  •  Etc. 

We  further  define  the  x(u,w)  in  the  rectangular  area  as: 
O^u^n,  G^,w^m,  and  n,m  is  the  positive  integer.  Further¬ 
more,  taking  the  ui  and  u^  values  in  formula  (9)  wherein  i=0,l, 
n,  j=0, 1,  * '  *  ,m,  then 

uw**  P(u)B,qT(u>)  (12) 


(10) 

(11) 
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In  the  formula: 


00 

00. 

!  01 

01.  . 

|  0m 

0m. 

00. 

O 

o 

* 

01. 

01..  . 

0  m« 

0m„ 

10 

10. 

11 

11.  . 

1  m 

10. 

10kW 

11. 

11..  . 

1m. 

1m.. 

.  !  .  . 

«0 

»0. 

1 

nm 

nm. 

nO. 

j 

i 

nmm 

i 

F,(ar), 

(ft 

-  ») 

F,(b,), 

(  l 

=  n 

Gt  (o/). 

(ft 

=  «  + 

i) 

Gt(b,), 

(  l 

=  j  + 

i) 

F  i  (<*/)» 

( ft 

-  i  + 

2) 

<7*(  “)  =  ' 

Ft(br), 

(  / 

=  i  + 

2) 

<?i(ai). 

(ft 

=»  i  + 

3) 

G,(b;), 

(1 

=  j  + 

3) 

0, 

m t 

(11 

0, 

nib 

(2) 

(13) 


(14) 


Key :  1 .  Other 
2 .  Other 


Moreover,  there  is: 


(  h  "  1*  2,—,Zn  +  2, 

<  i  -  2  (  u  )  +  1 , 

>  /  -  2 (w)+  1 , 


/  *  1,  2 ,  —,2m  +  2, 

a»=*  «-(»), 

b,=  ui  -(ui). 


Formula  (12)  is  then  the  divariant  interpolation  function  de¬ 
fined  on  the  rectangular  area  of  the  unit  square  mesh.  Each 
order  of  the  partial  derivative  in  the  Bx  matrix  is  very  dif¬ 
ficult  to  give.  Here,  we  use  the  smallest  second  power  method 
to  seek  B  [141  . 

X 

2.  Automatic  Mesh  Division 

The  form  and  nodal  point  coordinate  of  the  element  both 
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utilize  automatic  mesh  division  (this  only  requires  the  input 
of  a  small  amount  of  information) .  Its  process  is  divided  into 
two  steps:  first  is  to  determine  any  transverse  tangent  plane 
and  divide  it  into  tetrahedrons  according  to  any  type  of  stand¬ 
ard,  and  second  is  to  join  the  vertical  direction  to  make 
hexahedron  elements. 

The  cutting  up  of  the  transverse  tangent  plane  depends  on 
the  shape  and  position  of  the  cooling  holes;  for  a  solid  blade, 
we  can  cut  equally  or  according  to  proportion  along  the  arc 
length;  for  a  hollow  circular  holed  profile,  we  can  use  the 
center  of  the  circle  as  the  basic  standard  to  cut  along  the  cir 
cumference;  for  a  specially  shaped  holed  profile,  we  can  make 
a  curved  surface  fitting  similar  to  that  of  the  outer  form  be¬ 
cause  the  inner  walls  of  the  holes  also  have  curved  surfaces. 
Fig.  5.  contains  examples  of  various  types  of  transverse 
tangent  planes. 


Fig.  5. 


Cut  on  Profile's  Transverse  Tangent  Plane 

Key:  1.  Solid  blade 

2.  Hollow  blade 

3.  Circular  holes 

4.  Specially  shaped  holes 

VI.  Calculation  Results 

1.  Assessment  of  Examples 

(1)  Using  a  trapezoidal  sheet  with  varying  thickness  as  the 
example,  we  calculated  its  centrifugal  force  and  aerodynamic 
distribution  load;  using  a  hollow  cylinder  as  the  example,  we 
calculated  its  temperature  field  and  thermal  stress;  they  all 
had  analytical  solutions.  When  we  compared  it  with  the  calcula¬ 
tion  results  of  this  method  the  patterns  and  data  all  tallied 
and  satisfied  the  engineering  requirements. 

(2)  When  comparing  the  calculations  of  the  large  twisted 
blade  with  the  tests,  using  the  semi-solution  analytic  method  of 
reference  [5],  this  method's  error  for  the  large  twisted  (>20°) 
blade  was  relatively  large.  When  tested  at  n=17000  revolutions/ 
minute,  we  used  the  direct  measurement  method  yet  the  surface 
measured  point  was  not  the  same  as  the  principal  stress  dir¬ 
ection  and  thus  there  was  also  an  error. 

As  shown  in  fig.  6.,  the  calculation  results  of  this  paper 
basically  fall  within  the  tests  and  semi-analytic  method. 
Moreover,  this  method's  curve  is  smooth  and  is  therefore  con¬ 
sidered  believeable. 
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<  6 ) 

Fig.  6.  Test  and  Calculation  Curve  of  a  Large  Twisted  Blade 

Key:  (a)  Blade 

(b)  Data  on  measured  points 

1 .  Measured  point 

2 .  Method 

3 .  Test  point 

4.  Finite  element  method 

2 

5.  Kilograms/millimeter 

2.  Shrouded  Hollow  Air-Cooled  Blade 

(1)  Calculation  of  centrifugal  force 

For  the  calculated  objects  shown  in  fig.  1.,  in  order  to 
economize  on  the  internal  memory,  we  used  a  four  holed  tangent 
plane  to  carry  out  research.  For  the  stress  distribution  on  this 
blade's  .11  and  IV  tangent  planes  see  fig.  7.  The  high  stress  area 
is  on  the  tail  edge  of  the  tangent  plane  and  here,  besides  having 
a  centrifugal  force,  there  is  also  added  a  centrifugal  bending 
moment.  The  stress  distribution  of  the  2,3,4  points  on  each 
transverse  tangent  plane  along  the  blade's  high  speed  direction, 
as  shown  in  fig.  8.,  coincides  with  the  results  of  most 


calculation  methods.  The  stress  gradients  along  the  sides  of 
each  hole  are  not  large.  This  is  because  the  axis  line  of  the 
cooling  holes  are  basically  the  same  as  the  direction  of  the 
centrifugal  force  and  the  added  bending  moment  is  very  small. 
Naturally,  the  stress  value  near  the  small  holes  is  greater  than 
that  of  the  large  holes.  The  stress  value  on  the  shrouu  is  then 
very  small  and  the  shroud  is  equivalent  to  adding  a  concentrated 
mass  load  on  the  wingtip. 


Fig.  8.  Stress  Distribution  Along  Elade  Height 

Key:  1.  Tangent  plane  R 

2 

2.  (Kilograms/millimeter  ) 

3.  Center  of  gravity 

(2)  Aerodynamic  force 

The  blade's  stress  distribution  caused  by  aerodynamic  force 
is  basically  identical  to  the  stress  distribution  pattern  on  a 
varying  thickness  blade.  Yet,  in  an  actual  turbine  cascade,  the 
aerodynamic  force  changes  the  blade's  leading  and  trailing  edges 
as  well  as  the  blade  height  direction.  It  requires  the  use  of  a 
three  dimensional  flow  field  for  calculation  because  it  is  dif¬ 
ficult  to  accurately  give  the  aerodynamic  boundary  value.  This, 
then,  effects  the  accuracy  of  stress  calculations. 

(3)  Here,  we  will  only  give  the  temperature  field  for  the 
TT  and  f\T  tangent  planes  as  shown  in  fig.  9.  As  regards  the 
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temperature  distribution  of  each  cooling  hole  side,  the  cooling 
effect  of  the  large  holes  is  better  than  that  of  the  small  holes 
The  temperature  changes  along  the  periphery  of  each  hole  is 
not  large,  yet  the  normal  temperature  gradient  changes  on  the 
side  of  the  holes  is  large.  The  temperature  difference  between 
the  side  of  the  holes  and  the  profile's  outer  surface  is  about 
300°C.  The  high  temperature  area  is  on  the  leading  edge  of  the 
profile  and  the  tail  edge  location  because  the  leading  edge 
comes  in  direct  contact  with  the  gas  flow  and  the  tail  edge's 
dissipation  of  heat  is  difficult;  the  heart  shape  of  the  blade 
becomes  a  low  temperature  area ,  the  temperature  on  the  outer  sur¬ 
face  is  relatively  high,  the  temperature  field  along  the  height 
of  the  blade  is  related  to  the  temperature  and  flow  quantity  of 
the  air-cooled  gas  flow,  and  the  temperature  on  the  side  of  the 
blade  hole  is  close  to  the  temperature  of  the  air-cooling  gas 
flow.  The  curves  in  fig.  10.  are  the  temperature  changes  of 
points  1,3,2  from  the  blade  root  to  the  blade  tip.  Their  temper¬ 
ature  changes  are  all  about  250°C  and  the  shape  of  each  line  is 
similar;  the  temperature  on  the  blade  tip's  tail  edge  is  very 
high  and  this  seems  true  for  most  turbine  blades.  Cracks  and 
excessive  thermal  burning  are  often  produced  in  this  local 
position. 

(4)  Thermal  stress 

The  changes  of  the  thermal  stress  field  are  basically 
identical  to  the  change  patterns  of  the  temperature  differences. 
The  thermal  stress  on  each  tangent  plane  in  the  side  edge  of  the 
air-cooling  holes  and  outer  surface  is  very  high.  Although  the 
temperature  on  the  side  of  the  hole  is  not  high  yet  it  very 
easily  produces  cracks.  The  temperature  on  the  tail  edge  of  the 
blade  which  is  very  high  lowered  the  material's  permissable 
stress  value.  In  brief,  a  relatively  high  temperature  gradient 
area  causes  the  thermal  stress  gradient  to  also  be  relatively 
high. 


135 


Fig.  10.  Temperature  Distribution  Along  Blade  Height 

Key:  1.  Tangent  plane  R 
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A  STUDY  OF  THE  METHODS  FOR  MEASURING  RADAR'S  ECCM  PERFORMANCE 

Li  Nengjing 

(PLA  Air  Force  Laboratory) 

Abstract 

This  article  reviews  the  methods  for  measuring  radar’s 
ECCM  capabilities  and  points  out  a  set  of  commonly  used  measur¬ 
ing  formulas.  The  basic  parts  of  these  formulas  which  are  com¬ 
posed  from  the  radar's  principal  technical  parameters  repre¬ 
sent  the  radar's  potential  ECCM  capabilities.  The  supplementary 
parts  of  the  formulas  which  are  composed  from  the  indices  of 
the  radar's  various  ECCM  technical  measures  represent  the  quality 
level  of  these  ECCM  measures.  The  use  of  this  set  of  formulas 
can  calculate  the  expressed  numerical  values  of  the  radar's 
anti-passive  jamming  capabilities,  anti-active  jamming  capabil¬ 
ities  and  integrated  ECCM  capabilities.  This  paper  gives  ex¬ 
amples  and  discusses  the  application  of  the  calculation  results 
for  radar  systems  analysis  and  overall  design. 

I.  Preface 

Radar's  ECCM  capability  is  one  of  the  most  important  func¬ 
tions  of  radar  yet  up  to  the  present  there  is  still  not  a  com¬ 
monly  used  measuring  method. 

In  the  past,  there  were  two  commonly  used  means  of  express¬ 
ing  radar's  ECCM  capabilities.  One  was  the  listing  of  the  radar’s 
various  ECCM  technical  measures  and  their  performance  indices. 

For  example,  whether  the  radar  could  quickly  convert  its 
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frequency,  how  wide  is  its  frequency  conversion  band,  whether 
the  radar  has  moving  target  indication  (MTI)  and  how  much  is 
its  visible  coefficient  in  the  complex  waves  (SCV)  etc.  This 
type  of  expressed  formula  is  not  only  overelaborate  but  the 
ECCM  measure  items  and  index  itself  are  not  able  to  accurately 
measure  the  radar's  real  ECCM  capabilities.  Taking  the  radar's 
anti-passive  jamming  capabilities  as  an  example,  it  is  related 
to  the  SCV  value  as  well  as  the  radar's  resolution  or  pulse 
volume. 

The  second  type  of  formula  is  the  use  of  the  level  of  the 
radar's  check  or  measurement  performance  decline  under  specific 
interference  backgrounds  to  show  the  radar's  ECCM  capabilities. 
For  example,  under  a  specific  interference  power  spectrum 
density,  how  far  the  radar’s  defense  was,  how  much  the  radar's 
tracking  precision  of  a  certain  target  decreased,  etc.  The 
numerical  values  obtained  from  this  type  of  formula  which  are 
connected  to  the  interference  conditions  and  target  character¬ 
istics  cannot  independently  evaluate  the  radar's  ECCM  capabil¬ 
ities. 

J.  L.  Johnston  suggested  the  use  of  the  ECCM  improved 
factor  EIF  to  indicate  the  radar ' s  ECCM  performance  f 1 ] .  The 
definition  of  the  EIF  is 


In  the  formula,  (S/J)  is  the  radar's  output  signal  interference 

°  power  ratio  before  adopting  ECCM  measures; 

(S/J),  is  the  radar's  output  signal  interference 
*  power  ratio  after  adopting  ECCM  measures. 

This  formula  is  suitable  for  various  types  of  interference 
and  ECCM  measures  and  its  universality  is  relatively  strong.  Yet, 
it  is  only  suitable  for  measuring  the  above  type  of  radar  ECCM 
measures  on  the  integrated  ECCM  performances  of  several  types 


of  ECCM  measures  but  cannot  measure  the  entire  radar  system's 
ECCM  capabilities.  Because  of  certain  technical  parameters  of 
radar,  for  example  the  transmission  power  and  antenna  gain,  they 
are  obviously  not  able  to  belong  to  the  radar  's  ECCM  measures 
yet  they  determine  the  fundamental  factors  of  the  entire  radar's 
ECCM  capabilities. 

Because  of  this,  we  attempted  to  find  a  new  method  to  meas¬ 
ure  the  radar's  ECCM  capabilities.  The  demands  for  this  were: 
first,  to  only  rely  on  the  radar's  technical  parameters;  second, 
having  universality,  it  could  measure  the  radar's  capabilities 
to  counter  various  major  types  of  interference.  Based  on  these 
demands,  in  1978  for  the  first  time  I  proposed  an  expression 
formula  for  radar  integrated  ECCM  capabilities,  that  is: 

PT  B1GF.FC[2]  (2) 

O  1  AS 


In  the  formula 

P  is  the  radar's  mean  emission  power; 

Tq  is  the  radar's  observation  time  of  the  target; 

B.  is  the  radar  system  's  instantaneous  band 
width; 

G  is  the  radar  antenna's  gain; 

Fa  is  the  radar  antenna's  quality  factor; 

F  is  the  radar  signal's  quality  factor. 

This  paper  will  further  prove  and  investigate  this  method. 

II.  An  Analysis  of  the  Interference  Type  and  the  Fundamental 
Measures  of  Radar's  ECCM  to  Bring  Forth  an  Expression 
Method  for  Radar's  ECCM  Capabilities 

In  order  to  bring  forth  an  expression  method  for  radar's 
ECCM  capabilities,  we  first  observed  the  radar's  fundamental 
ECCM  measures.  See  table  1. 
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Table  1 . 


Type  of  Interference 


Radar  1 s  Fundamental  ECCM  Measures 


1. 


Emission  type  interference 
(active  jamming) 


(a) 


Increase  the  radar's  signal 
power . 


(1)  Inhibition  such  as 
noise  interference 
or  other  modulation 
wave  form  interfer¬ 
ence. 


(b)  Cause  a  difference  between 
the  radar  1 s  signal  carrier 
frequency  and  interference  - 
resolved  from  the  frequency. 

(c)  Raise  the  radar  antenna's 
orientation  -  resolved  from 
the  space  direction  angle. 


(d)  Have  the  generally  best  filter 
so  as  to  find  the  greatest 
output  signal  interference 
power.  Use  a  matching  filter 
for  white  noise  interference; 
then  join  the  most  effective 
filter  of  this  type  of  in¬ 
terference  wave  form  for  the 
non-white  noise  interference  - 
resolved  from  the  signal  wave 
form. 


(2)  Deceiving  such  as 
response  jamming. 


2.  Reflection  interference 
(passive  jamming)  such 
as  interference  tinsel 
cord  and  non-man  made 
surface  features,  clouds 
and  rain,  and  ocean  waves. 


(a)  Design  a  complex  radar  signal 
so  that  the  jammer  finds  it 
difficult  to  duplicate  - 
resolved  from  signal  wave. 

(b)  Design  a  radar  operations 
system  so  that  the  inter¬ 
relation  between  the  radar 
signal  wave  and  space  direction 
angle  are  smallest  (dealing 
with  angle  deception  response 
interference) . 

(a)  Use  different  speeds  for  the 
radar  target  and  interference 
complex  wave  (the  Doppler  fre¬ 
quency)  ,  signal  selected  from 
complex  waves  -  resolved  from 
frequency. 


i 

f 

i 

i 
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Type  of  Interference 


Radar's  Fundamental  ECCM  Measures 


3.  Absorbing  interference 
such  as  electromagnetic 
absorption  coating. 


(b)  Raise  the  radar  antenna 

orientation  and  reduce  the 
radar  signal  time  width  so 
as  to  decrease  the  strength 
of  the  received  complex  wave, 
-  that  is,  raising  the 
resolution  for  the  space 
direction  angle  and  time. 

(a)  Enlarge  the  radar  signal 
power . 


From  an  analysis  of  the  above  brief  table,  besides  the  spec¬ 
ial  demands  for  dealing  with  angle  deception  response  interfer¬ 
ence  with  signal  wave  form  and  space  direction  angle,  there  are 
two  basic  radar  measures  for  ECCM:  one  is  power  and  the  other 
is  resolution.  Because  of  this,  it  can  be  inferred  that  by  using 
the  radar's  power  parameter  and  resolution  parameter,  we  can 
indicate  the  radar's  fundamental  ECCM  capabilities. 

To  represent  the  radar  power's  technical  parameters,  we 
should  use  the  radar's  mean  transmission  power  P. 

The  radar's  resolutions,  for  example  those  mentioned  in 
table  1,  have  the  aspects  of  frequency,  space  directional 
angle,  time  and  signal  wave  form.  From  the  signal  analysis 
theory  it  can  be  known  that,  in  essence,  the  wave  forms  resol¬ 
ution  is  the  resolve  of  the  signal  in  the  time  domain  and  fre¬ 
quency  domain.  Because  of  this,  the  radar's  resolution  is 
mainly  concerned  with  four  areas  of  overall  resolution:  the 
time  domain,  frequency  domain,  space  directional  angle  and 
angle  of  elevation. 

Without  the  simple  parameters,  we  can  still  express  the 
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radar's  integrated  resolution.  Yet,  we  know  that  the  radar  four 
dimension  blur  function  brought  forth  by  H.  Urkowitz  [3]  can 
describe,  on  a  relatively  wide  foundation,  the  integrated  resol¬ 
ution  of  radar  in  space,  time  domain  and  frequency  domain. 

The  function  is: 


A  (.It,  ft,  -/,)F*(9,  4>)Ft(0  -0.) 

-  oo 

X  (4>  -b<)'le-il’,‘dfdQdb  (3) 


In  the  formula, 

U (f )  is  the  radar  signal’s  envelope  complex  frequency 
spectrum; 

F(0,£)  is  the  radar  antenna's  direction  fiqure; 
t^  and  f^  are  the  time  difference  and  frequency  difference; 

$=sin  a  is  the  sine  value  of  directional  angle  a; 

<b-sin  £  is  the  sine  value  of  elevation  angle  £  ; 

6,  and  are  the  difference  of  &  and  the  difference  of  <t. 
a  d 

When  the  antenna  band  width  is  far  greater  than  the  signal 
band  width  so  that  there  is  no  apparent  interaction  between  the 
antenna  and  signal. 


A  (ft,  ft,  0*.  4><)  *  A  (ft,  ft)  A  (9 1, 


(4) 


This  is  the  product  of  the  four  dimensional  blur  function  which 
can  be  expressed  as  the  signal  blur  function  and  the  antenna 
direction  blur  function.  Formula  (4)  is  applicable  for  most 
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i 


radar . 

The  signal  blur  function  ACt^.f^)  is: 


*  p Um{f)UU 

”f  “(*>«*(<  +(i)eit‘lddt  (  5  1 


In  the  formula,  u(t)  is  the  signal's  envelope  time  function. 

The  antenna  direction  blur  function  A(8,,  ^Jis: 

a  a 

oo 

- JJf-o,  4>)FC(e-fl4),(*-WV9«/*  (e) 


This  function  is  the  self  correlation  function  of  the  antenna 
direction  figure.  Because  most  of  the  radar  antenna  direction 
figure  is  close  to  the  sinc  function,  the  self  correlation  func¬ 
tion  of  the  sinc  function  is  still  the  sin  function.  Because  of 
c  c 

this,  the  antenna  direction  figure  is  clore  to  its  self  correl¬ 
ation  function.  We  can  use  the  antenna  direction  figure  to  sub¬ 
stitute  for  the  antenna  direction  blur  function. 

In  order  for  the  resolution  to  be  able  to  be  expressed  by 
a  parameter,  we  used  the  main  peak  of  blur  function  modular 
value  IaI  when  it  decreased  to  a  specified  range.  For  example, 
when  it  decreased  to  the  main  peak  and  the  highest  value  (after 
normalizing  this  value  was  1)  was  -6  decibels,  the  width  of 
the  main  peak  (when  the  blur  function  is  one  dimensional) ,  or 
the  sectional  area  (when  the  blur  function  was  two  dimensional) , 
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or  sectioned  volume  (when  the  blur  function  is  even  more  multi¬ 
dimensional)  acts  as  the  parameter.  This  parameter  can  be  called 
the  resolution  unit. 

We  call  the  sectional  area  of  the  two  dimensional  blur 
function  lACt^f^)]  whose  main  peak  decreased  6  decibels  the  radar 

signals  distance  and  frequency  integrated  resolution  unit.  The 
sectional  area  of  the  antenna  direction  figure  )  F  ( 0 ,  <#>)  |  whose 
peak  decreased  6  decibels  is  called  the  radar  antenna's  space 
angle  resolution  unit. 

It  is  inferred  from  this  that  the  four  dimensional  volume 
obtained  when  the  four  dimensional  blur  function 

main  peak  decreased  6  decibels  can  be  defined 
as  the  radar's  integrated  resolution  unit.  When  R^  is  used  to 

represent  this  value,  the  larger  the  R^,,  the  more  error  shown 
for  the  radar's  integrated  resolution.  Therefore,  the  precise 
description  of  the  radar  resolution  parameter  is  — . 

i 

Radar's  basic  measure  for  anti-passive  jamming  is  raising 
the  radar's  resolution.  Therefore,  the  description  of  the  para¬ 
meter  of  the  radar's  integrated  resolution  (^)  can  be  used  to 

(RT) 

express  the  radar's  basic  anti-jamming  capabilities. 

Padar ' s  basic  methods  for  anti-inhibitory  active  jamming 
are  to  increase  the  energy  of  the  target  echo  as  much  as  possible, 
decrease  the  interference  power  coming  in  from  the  antenna  as 
much  as  possible  and  filter  out  the  interference  that  has  already 
been  transmitted  in  so  that  the  power  ratio  of  the  output  signal 
and  interference  will  be  maximal.  To  realize  these,  we  must 
rely  on  radar  with  great  emission  power  and  at  the  same  time 
depend  on  radar  with  high  resolution.  If  the  latter's  emission 
power  can  be  concentrated  on  the  target,  then  both  the 
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interference  input  and  output  will  decrease.  Therefore,  the 
product  of  the  power  parameter  and  resolution  parameter  which  is 
(P/R,p)  can  be  used  to  express  the  radar's  basic  anti-actjve 
jamming  capabilities. 

As  regards  the  response  type  interference,  besides  havinq 
specially  interrelated  demands,  we  can  also  depend  on  the  radar's 
high  resolution.  As  regards  absorption  interference,  it  is  also 
necessary  for  the  emission  power  to  be  concentrated  on  the 
target.  This  is  identical  to  dealing  with  active  jamming.  From 
a  wide  general  perspective,  the  product  of  the  power  and  resol¬ 
ution  parameter  (P/RT)  can  be  used  to  express  the  radar's  funda¬ 
mental  integrated  ECCM  capabilities  in  dealing  with  various  types 
of  interference. 

Yet,  the  calculation  of  the  R^  numerical  value  is  complex 
and  an  inferred  simplified  formula  should  be  substituted. 

Firstly,  we  established  formula  (4)  for  most  radar  and 
because  of  this, 

I A  0*  /*  0*  40I  =  U  (0,  /«)!l  A  (9*  <M  (■  7  > 

In  this  way,  RT  is  then  equal  to  the  four  dimensional  volume 
formed  from  the  3  decibel  modular  value  sectional  surface  of  the 
signal's  blur  function  and  the  3  decibel  modular  value  sectional 
surface  of  the  antenna  direction  blur  function. 

There  are  three  types  of  radar  signal  blur  function  modular 
value  forms  (abbreviated  as  the  three  dimensional  blur  figure) : 
nail  board  form,  thumbtack  form  and  knife-edge  form  [4-5], 

When  the  radar  signal  is  a  pulse  train,  the  three  dimensional 
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Fig.  2.  Plane  Blur  Figure  of  Linear  Frequency  Modulation  Pulse 
Train 

We  know  from  figures  1  and  2  that  the  approximate  value  of 
this  type  of  signal  blur  function  modular  value  main  peak  with 
a  3  decibel  sectional  area  is: 

*  J  _L_  J*_  1 

4  NT  Bs  “  4 

When  the  radar  signal  is  formed  by  a  single  relatively  long 
pseudohatted  code,  then  the  three  dimensional  blur  function  of 
this  type  of  signal  is  a  thumbtack  form.  Its  plane  blur  figure 
is  shown  in  fig.  3. 


Fig.  3.  Plane  Blur  Figure  of  a  Single  Pseudohatted  Code  Signal? 
"t*  is  the  Total  Length  of  the  Signal 

From  fig.  3  we  know  that  the  approximate  value  of  this  type 
of  signal  blur  function  modular  value  main  peak  with  a  3  decibel 
sectional  area  is: 


_3_  I  3  *  

4  4  l  qUs 

Given  that  at  this  time  the  signal's  sustained  time  is  T  ,  then 

o 

the  total  length  of  the  code  is  /T"o. 

When  the  radar  signal  is  a  single  linear  frequency  modula¬ 
tion  pulse,  the  three  dimensional  blur  figure  of  the  signal  is 
a  knife-edge  form.  See  fig.  4  for  its  plane  blur  figure. 


Fig.  4  Plane  Blur  Figure  of  a  Single  Linear  Frequency  Modulation 
Pulse 

This  type  of  signal  can  only  accurately  determine  the  com¬ 
bined  value  of  an  unknown  target's  distance  and  speed.  This  com¬ 
bined  value  is  the  oblique  axis  in  the  plane  blur  figure.  Yet, 
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we  cannot  precisely  know  what  the  distance  and  speed  are.  There¬ 
fore,  at  the  least,  during  actual  application  we  must  transmit  a 
pair  of  linear  frequency  modulation  signals.  This  pair  of 
signals  have  opposite  linear  frequency  modulation  slopes  as 
shown  in  fig.  5. 


Fig.  5.  Plane  Blur  Figure  of  a  Pair  of  Linear  Modulation 
Frequency  Pulses 


At  this  time,  a  target  will  produce  two  echos  on  the  indicator? 
the  intermediate  value  of  the  two  echos  is  the  accurate  dis¬ 
tance  and  the  interval  of  the  two  echos  represents  the  target ' s 
speed.  The  equivalent  plane  blur  figure  of  the  radar  signal 
composed  by  this  pair  of  linear  frequency  modulation  pulses  is 
the  overlapping  part  of  the  two  directional  slope  ellipse  which 
is  the  oblique  line  part  in  fig.  5.  This  area  is  approximately 
equal  to  ^x(  r  •~)  ..  When  considering  T  =2  T  ,  then  this  area  is 

to  O 

i  s 

equal  to  *y~  • 
o  s 

We  know  from  the  above  analysis  that  the  value  of  the  com¬ 
monly  used  radar  signal  blur  function  modular  value  main  peak 
with  a  3  decibel  sectional  area  is: 


( * ) 
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In  the  formula,  k^  is  the  coefficient  and  its  approximate  value 
is  3/4  ~1. 

The  antenna  direction  blur  function  can  use  the  antenna 
direction  figure  for  an  approximation.  Therefore,  the  approximate 
value  of  the  antenna  direction  blur  function  modular  value  main 
peak  with  a  3  decibel  sectional  area  is: 


4  a°e»  (9) 


In  the  formula,  a  and  l  are  the  3  decibel  widths  of  the  antenna 

o  o 

direction  figure's  directional  angle  and  elevation  angle. 


Thus,  we  can  now  obtain  the  approximate  value  of  RT: 


Rr**(k*  ^  4  a°&°  )  — 


4 


aoe«_ 

TtBs 


(10) 


Because  the  antenna  gain  G  can  be  approximately  expressed  as 


G  = 


a,e0 


(11) 


Therefore 


l 

T,B,G 


(12) 


Because  the  constant  coefficient  is  of  no  importance  it  can 
be  eliminated  in  normalization  and  therefore  the  radar's  inte¬ 
grated  resolution  parameter  (■=—■)  can  be  simplified  and  expressed 

T 

as  T  B  G. 
o  s 
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p 

Using  (Rt) ,  the  expressed  radar  fundamental  integrated  ECCM 
capabilities  can  be  simplified  as  : 


PT,B,G 


(13) 


III.  The  Expression  of  Radar  Fundamental  ECCM  Capabilities 
Derived  From  the  Signal ' s  Interference  Power 

The  expression  of  radar's  fundamental  ECCM  capabilities  ob¬ 
tained  in  the  last  section  can  also  use  another  method  which  is 
derived  from  the  radar's  transmitted  signal  interference  power 
ratio. 

If  (^) j  represents  the  power  ratio  of  a  radar  target  signal 
and  passive  jamming  in  a  radar  input  area  and  represents  the 

power  ratio  of  *-he  radar's  output  signal  and  passive  jamming,  then 


In  the  formula 


d  is  the  target's  radar  sectional  area? 

dQ  is  the  radar  sectional  area  of  the  passive 
jamming. 


o. «  |  ct/f’gn,  (I-*) 


In  the  formula, 

c  is  the  electric  wave  propagation  speed; 

R  is  the  distance  from  the  passive  jamming 

area  center  to  the  radar  station; 

Q  is  the  three  dimehsional  angle  value  of 
the  radar  antenna  direction  figure; 

*1  is  the  reflex  coefficient  of  the  passive 
jamming. 
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Because 


T 


Therefore 


/  S  N  o,n,G_  (in) 

\  C  /(  2 xcVjV 


and 


m-'(sA 

i 

In  the  formula,  I  is  the  modified  factor. 

If  the  speed  distribution  of  the  interference  object  is 

basically  uniform,  its  Doppler  frequency  band  width  is  B  ,  its 

c 

signal  frequency  band  width  is  B  ,  B  ^  B  and  the  two  are 

s  s  c 

B 

overlapping,  then  the  theoretical  value  of  I  is  _c  .  In 

B 

s 

real  situations,  if  the  interfering  frequency  band  and  signal  band 
are  possibly  not  overlapping  or  not  completely  overlapping  then 
we  can  define  Bc  as  the  equivalent  interference  band  width  which 
is: 
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(18) 


B.=IBt 


We  already  know  that 


(19) 


and  therefore 


I 


**  -r.t% 


(20) 


Because  of  this 


(  S  \  _ 

\  C  )*  2*c\R2' 


kx(T,B,G) 


(21) 


In  the  formula. 


.  _ o.B. 

1  2nc\R2 


(22) 


It  is  composed  by  the  constant  coefficient  2  If  c  and  target 

2 

paramenter  ^  as  well  as  the  interference  conditions  and  Rc . 

When  deciding  the  size  of  the  radar  output  signal-interfer¬ 
ence  ratio  in  passive  jamming,  we  also  decided  that  the  radar’s 
technical  parameter  for  the  anti-passive  jamming  capabilities 
size  is  TqBsG.  This  is  a  simplification  of  the  radar  integrated 

resolution  parameter  -4— derived  in  the  last  section. 
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Let  us  take  a  further  look  at  the  situation  of  anti-inhib¬ 
ition  active  jamming. 

Letting  (j)Q  represent  the  output  signal  interference  power 
ratio  when  the  radar  receives  active  jamming  and  if  the  inter¬ 
ference  power  spectrum  is  basically  uniform  (this  is  a  relatively 
common  interference  condition) ,  then  the  radar  selected  and 
signal  matched  receiver  at  this  time  is: 


In  the  formula, 

Eg  is  the  energy  of  the  target's  echo  signal; 
P_  is  the  received  interference  power. 


Therefore 


-  PT'G'X'-o, 
ts~  64ns/?i 


(  > 


71  is  the  radar's  wavelength; 

R  is  the  distance  between  the  target  and  the 
radar  station. 


In  the  formula. 


In  the  formula. 


Pj  is  the  jammer's  reflection  power; 

Gt  is  the  jammer's  antenna  gain; 
u 

A  (9)  is  the  radar  antenna's  equivalent  receiv- 
~  ing  area  for  the  interference  direction; 

Yj  is  the  polarization  cc  .^ficient; 

R_  is  the  distance  between  the  jammer  and  the 
radar . 


If  the  jammer's  location  direction  is  basically  the  same  as 
the  target,  then 


A,(  0  )  = 


£X_* 
4  K 


(26) 


After  substituting  in  formula  (25) : 


p  PiGtYi’k'G  _ 


(27) 


When  formulas  (24}and  (27)  are  substituted  into  formula  (23): 


(  S\ _ 

7  °  4si(^j")P/G/Y/ 


=  k.(PT  aBsG) 


(28) 
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In  the  formula. 


k.  = 


■It) 


P/GjYj 


(29) 


It  is  composed  of  target  parameters  and  R  ,  and  inter¬ 
ference  conditions  P_,G_,  YT  and  RT. 

J  U  J  J 

Therefore,  when  deciding  the  size  of  the  radar  output  signal 
interference  ratio  under  active  jamming,  we  also  decided  that  the 
radar's  technical  parameter  for  anti-active  jamming  capabilities 
size  is  PTQBgG.  This  is  a  simplification  of  the  radar  power  and 
resolution  parameter  product  )  derived  in  the  previous 

r< 

T 

section. 

Now  we  will  take  a  further  look  at  the  radar  output  signal- 
interference  ratio  condition  when  two  types  of  interference  exist 
simultaneously.  Let  represent  the  output  signal -jamming 

ratio . 

_ _ l _ 

1  +1_  1  .  - 
kJ.BsG  ktPT„BsG 

( PT,BSG )  (30) 

If 


1 


(31) 
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Then 


(,s+cb^PT-B’G> 


(32) 


which  is  the  same  as  when  there  is  only  active  jamming. 

When  we  employed  a  representative  technical  parameter  to 
calculate  the  value  of  (^2^  ,  we  could  prove  that  formula  (32) 

V 

is  established  in  a  rather  wide  range.  Therefore,  formula 
PTQBgG  not  only  expressed  the  radar's  fundamental  anti-active 

jamming  capabilities  but  it  was  also  able  to  express  the  radar's 
integrated  fundamental  ECCM  capabilities  for  anti-passive  jam¬ 
ming  and  anti-active  jamming. 

IV.  Supplementary  Factors  for  the  Fundamental  Expression 

The  already  derived  radar  fundamental  ECCM  capabilities  ex¬ 
pression  measures  the  radar's  fundamental  or  potential  ECCM 
capabilities.  If  we  want  to  measure  the  radar's  ECCM  capabilities 
more  completely  it  is  also  necessary  to  have  the  usefulness  or 
quality  level  of  the  radar's  various  ECCM  measures  act  as 
supplements.  Below  we  will  mention  a  group  of  supplementary 
factors . 

1 .  Quick  Changing  Carrier  Frequency 

Because  the  radar's  integrated  four  dimensional  blur  function 
only  takes  into  consideration  the  signal's  complex  envelope  and 
has  not  yet  considered  the  radar's  carrier  frequency,  it  only 
expresses  the  resolution  of  the  signal  carrier  frequency  when  in 
a  static  state.  The  changes  of  the  carrier  frequency  are  often 
used  in  the  radar  ECCM  measures  to  resolve  the  signal  and  active 
jamming.  Yet,  these  changes  must  be  quick  changes  within  the 
pulse  or  within  a  small  group  of  pulses  to  be  able  to  shake  off 
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an  enemy  location's  track  jamming  and  force  the  enemy  to  use 
wide  frequency  band  jamming.  Because  the  total  power  of  an 
enemy’s  jammer  is  fixed  or  limited,  the  jamming  power  spectrum 
enlarges  times,  the  jamming  power  density  then  drops  kj 
times  and  the  radar’s  output  signal  jamming  power  ratio  then 
rises  kj  times.  Because  of  this,  the  quick  change  measures  of 
the  carrier  frequency  are  equal  to  raising  the  radar's  trans¬ 
mission  power  to  k^  times. 

The  radar  allows  the  maximum  range  of  its  carrier  frequency 
quick  changes  to  be  the  radar  system's  instantaneous  band  width 
B^ .  Naturally,  radar  designers  have  made  every  effort  to  reach: 


Therefore,  after  the  radar  has  quick  changing  frequency 
measures,  the  expression  of  the  radar's  anti-active  jamming  cap¬ 
abilities  should  add  a  supplementary  k^  factor. 

2.  Side  Lobes 

When  we  defined  the  radar's  integrated  resolution  parameters 
we  only  considered  the  sectional  volume  when  the  integrated  blur 
function  main  peak  decreased  6  decibels.  Because  of  this,  the 
corresponding  signal  blur  function  and  antenna  direction  figure 
also  only  considered  its  main  peak  3  decibel  area.  Yet,  there 
are  many  side  peaks  around  the  main  peak  and  they  influence  the 
ECCM  capabilities  within  a  fixed  range.  This  will  be  discussed 
below. 

When  the  present  microwave  radar  antenna  direction  figure 
main  peak  (commonly  called  the  antenna  beam  main  lobe)  is  very 
sharp-pointed,  the  active  jamming  carried  out  from  the  side 
peak  (commonly  called  side  lobe  or  secondary  lobe)  is  already 
the  principal  jamming  method  of  the  radar  against  the  enemy. 
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The  side  lobe  level  has  already  become  a  technical  index  for 
measuring  the  radar's  ECCM  capabilities.  Therefore,  in  the  ex¬ 
pression  of  our  measuring  the  radar's  anti-active  jamming  cap¬ 
abilities,  it  is  necessary  to  introduce  the  supplement  factory 
Gg  to  represent  the  radar  antenna's  side  lobe  level. 


We  defined  factor  G  as: 

s 


(decibel) =101g 


(antenna  power  direction  figure 

(main  lobe  peak  value 
/ 


^antenna  power  direction  figure  \ 
(maximum  secondary  lobe  peak  value) 


) -25  (34) 


Defining  it  in  this  way  takes  into  consideration  the  present  radar 

technical  level.  The  antenna's  main  side  lobe  power  ratio  of  25 

decibels  can  be  viewed  as  a  medium  level  representative  value. 

When  it  is  higher  than  25  decibels  this  represents  the  antenna's 

high  grade  design  which  causes  the  G  to  have  a  positive  value. 

s 

On  the  contrary,  when  it  is  lower  than  25  decibels  this  represents 

a  poor  design  causing  the  Gg  to  have  a  negative  value.  This  type 

of  method  is  used  to  cause  the  G  factor  to  influence  the  numer- 

s 

ical  value  of  the  original  anti-active  jamming  capabilities 
fundamental  expression. 

The  G  only  has  one  value  for  the  needle-shaped  beam.  The 
Gs  should  include  two  factors  for  the  fan-shaped  beam,  that  is, 
a  horizontal  side  lobe  factor  GSH  and  a  vertical  side  lobe  factor 
G  .  These  two  factors  are  separately  added  after  calculating 


formula  (34) . 


G  =G  +G 
s  SH  sv 


When  the  radar  has  an  offsetting  device,  the  offsetting 

benefits  should  be  added  into  G  .  Sometimes  when  the  radar 

s 

signal's  blur  function  is  in  a  relatively  distant  area  of  its 


main  peak  there  exist  many  relatively  high  branch  peaks  (also  ■ 
called  grating  lobes)  and  side  peaks  (also  called  side  lobes) 
near  the  main  lobe. 


The  branch  lobes  or  grating  lobes  usually  lie  in  the  time 
domain  of  the  radar's  sustained  interference  and  outside  of  the 
frequency  domain.  Here,  we  will  put  off  discussion  on  this 
temporarily. 

The  signal  side  lobe  is  a  signal  processing  circuit 
which  must  be  inhibited.  Its  disadvantage  lies  in  a  weak  signal 
producing  interference  near  a  strong  signal.  Yet,  it  is  differ¬ 
ent  from  the  antenna  beam  side  lobe  in  that  it  does  not  belong 
to  a  hostile  interference  entrance.  The  demands  for  the  inhib¬ 
ited  signal  side  lobes  also  do  not  resemble  the  antenna  side 
lobes  direct  columns  which  are  the  radar ' s  ECCM  performance. 
Because  of  this,  we  do  not  take  the  signal  side  lobe  inhibi¬ 
tion  as  an  ECCM  factor  for  calculations. 

3.  Signal  Processing 

The  radar  resolution  described  by  the  radar  blur  function 
indicates  the  potential  abilities  of  the  radar's  resolve  signal 
and  if  the  potential  abilities  are  realized  it  is  still  necessary 
to  adopt  measures  in  signal  processing.  Firstly,  the  radar  re¬ 
ceiver  should  be  a  matching  wave  filter  because  the  radar 
signal's  blur  function  is  essentially  formed  from  a  group  of 
different  Doppler  frequency  target  echos  which  pass  through  the 
output  wave  combination  of  a  matching  wave  filter  within  the 
same  period  of  time  (4.5].  Yet,  the  radar  does  not  necessarily 
use  an  integrated  matching  wave  filter.  For  example,  generally 
in  dealing  with  passive  jamming  the  radar  does  not  use  a 
Doppler  frequency  shunt  circuit  matching  wave  filter  but  after 
it  is  in  a  single  circuit  matching  wave  filter  it  connects  the 
MTI  circuit.  Further,  in  order  to  deal  with  non-white  noise 
type  active  jamming  it  is  also  necessary  to  connect  special 
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filters.  The  quality  and  anti- jamming  usefulness  of  these  signal 
processing  circuits  are  supplementary  factors  which  we  must  con¬ 
sider. 


Firstly,  we  will  consider  the  quality  of  the  MTI.  B.D. 
Steinberg  has  already  defined  quality  factor  I  for  measuring  the 
quality  of  the  MTI ' s  technical  measures  [6], 

MTI  output  signal  stray 
j  _  wave  power  ratio 

MTI  input  signal  stray 
wave  power  ratio 

SCV  is  often  used  in  a  radar  technical  performance  expression 
to  substitute  for  I.  The  value  of  SCV  is  6  decibels  lower  than  I. 

In  the  anti-passive  jamming  capabilities  expression,  we 
defined  a  moving  target  to  indicate  quality  factor  PM: 

PM (decibels) =SCV-25  (37) 


Taking  the  constant  25  decibels,  we  first  consider  that  it 
is  able  to  represent  the  canonical  level  of  the  modern  radar  MTI 
visibility  coefficient  SCV  and  at  the  same  time  cause  it  to  be 
convenient  for  memory  and  use  with  the  previous  side  lobe 
factor  constant. 

Second,  the  signal  processing  factor  which  must  be  considered 
is  a  constant  false  alarm  processing  quality  factor. 

The  derivation  of  the  radar's  integrated  resolution  formula 
takes  the  radar’s  receiving  system  linearity  as  a  basis.  If  the 
strength  of  the  signal,  interference  or  both  surpasses  the  dynamic 
range  of  the  radar's  receiving  system,  this  causes  the  receiving 
system  to  be  saturated  or  obstructed  and  the  radar ' s  resolution 

166 


to  decrease  drastically.  Because  of  this,  the  constant  false 
alarm  is  a  necessary  radar  measure  for  various  types  of  anti- 
j  amming . 


Therefore,  the  action  of  the  constant  false  alarm  process¬ 
ing  in  the  anti- jamming  is  anti-saturation  and  after  constant 
false  alarm  processing  the  signal  interference  ratio  has  a  loss. 
Because  of  this,  we  defined  the  constant  false  alarm  processing 
supplementary  factor  in  the  radar's  ECCM  capabilities  expression 
as : 

U)  / \M\ 

Pf(#fl)“10!g(^J-25  (38) 

Key:  1.  (Decibels) 

In  the  formula, 

A  M  represents  the  allowable  increased 

multiple  in  the  radar  receiving  system's 
dynamic  range  after  adding  in  the  con¬ 
stant  false  alarm  processing  device: 

Lcf  represents  the  constant  fale  alarm  loss: 

The  selected  theorem  of  constant  25  decibels 
is  from  the  same  previous  formula. 


Noise  frequency  modulation  jamming  is  the  often  encountered 
wide  band  noise  type  active  jamming.  In  order  to  deal  with  this 
type  of  interference,  modern  radar  has  already  made  wide  use  of 
Dicke-Fix  (i.e.  wide-limited-narrow  circuit) .  We  defined  the 
Dicke-Fix  supplementary  factors  for  the  anti-active  jamming 
capabilities  as: 

8  (39) 

Key:  1.  (Decibels) 

4.  Antenna  Polarization 

The  radar's  comprehensive  resolution  discussed  previously 
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did  not  include  the  influence  of  antenna  polarization.  There¬ 
fore,  it  is  necessary  to  further  analyze  the  polarized  resol¬ 
ution.  Theoretically  speaking,  because  the  radar  target's 
scattered  echo  and  outside  interference  are  completely  independ¬ 
ent,  the  polarization  states  can  never  be  the  same.  Because  of 
this,  they  can  be  used  to  cause  the  antenna  to,  as  best  as  pos¬ 
sible,  have  polarization  matching  for  the  useful  echo  signals 
and  exert  the  mismatched  method  for  the  interference  to  obtain 
ECCM  results.  Yet,  in  reality,  it  has  good  results  for  rain  and 
cloud  interference. 

Because  raindrops  are  nearly  circular,  most  of  its  reverse 
scattering  for  the  circular  polarized  waves  are  reverse  rota¬ 
tional  circular  polarized  waves  but  aircraft  and  other  complex 
man-made  targets  are  not  like  this.  Therefore,  the  use  of  a 
circular  polarized  antenna  or  an  antenna  whose  polarization  can 
be  randomly  adjusted  can  benefit  a  signal -interference  ratio  of 
10-20  decibels.  Yet,  complex  structured  natural  objects  which 
produce  interfering  complex  waves,  for  example  surface  object 
complex  waves,  are  not  able  to  obtain  this  benefit. 

As  regards  man-made  interference,  if  active  jamming  comes 
from  a  jammer  or  passive  jamming  object  which  possesses  rela¬ 
tively  uniform  scattered  polarization  characteristics  (for  ex¬ 
ample,  interference  tinsel  cords  hanging  down  all  assume  a  hor¬ 
izontal  state) ,  then  the  use  of  a  polarized  adjustable  antenna 
or  the  polarization  offset  method  can  have  very  good  ECCM  re¬ 
sults.  Yet,  in  an  electronic  war,  the  enemy  actually  takes  pre¬ 
cautions  against  these  measures.  Its  jammer  is  usually  a  pair 
of  simultaneously  started  machines  which  separately  transmit 
two  types  of  opposite  polarized  waves.  Its  put  in  interference 
wire  is  also  designed  to  hang  down  with  one  half  horizontal 
and  the  other  half  vertical.  Because  of  this,  under  actual 
countermeasure  conditions,  the  radar  antenna  feed  polarization 
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which  can  be  adjusted  can  only  benefit  3  decibels.  If  the  radar 

does  not  have  this  changeable  polarization  measure,  the  W  value 

s 

is  0. 


5.  Dealing  with  Response  Type-Interference 

There  are  several  different  types  of  situations  in  which 
response  type  interference  disturbs  radar. 

Repetitious  entering  from  the  secondary  lobes  is  usually 
used  for  search  radar.  This  causes  response  pulses,  disturbance 
and  useful  hidden  target  signals  on  many  locations  and  at  dif¬ 
ferent  distances  on  the  indicator  frame.  The  measures  for  deal¬ 
ing  with  this  type  of  interference  are  mainly  inhibiting  the 
antenna  side  lobes  and  causing  the  radar  to  have  repeated  fre¬ 
quency  chatter. 

For  tracking  radar,  the  response  type  interference  mainly 
has  gate  pulling  and  counter -modulation  angle  deception.  The 
counter  measures  for  the  former  are  wave  form  differentiation 
and  repeated  frequency  chatter.  The  basic  measures  for  the 
latter  are  using  a  signal  wave  form  and  space  direction  angle 
non-interrelated  operating  system;  for  example,  the  single 
pulse  goniometer  system. 

We  have  already  discussed  wave  form  resolve  and  antenna 
side  lobe  isotopes  previously. 

The  repetition  frequency  chatter  measure  is  already  being 
widely  used  in  modern  radar.  Besides  counter  response  type  in¬ 
terference,  it  is  also  able  to  inhibit  a  nearby  station's 
synchronic  interference  and  is  also  able  to  overcome  the  blind 
speed  of  the  MTI  technique.  As  regards  the  signal’s  three  dimen¬ 
sional  blur  figure,  its  action  causes  the  branch  peak  of  the 
nail  board  form  blur  function  to  lower  and  the  lowered  multiple 
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is  equal  to  the  numeral  J  of  the  chatter  changes.  Shown  on  the 
radar's  reception  terminal,  it  also  surpasses  the  periodic 
response  signal's  lowered  output  power  J  times.  For  this,  we 
defined  the  repetition  frequency  chatter  measure’s  supplementary 
factor  for  the  radar's  ECCM  capabilities  as: 

(1) 

P/(#5l)-101g/  -  8  (40) 

Key:  1.  (Decibels) 

By  taking  a  constant  8  decibels,  because  the  modern  radar  J 
is  commonly  5  to  7,  8  decibels  represents  all  levels.  At  the  same 
time,  it  is  also  identical  to  formula  (39)  and  is  thus  easy  to 
remember . 

The  reason  that  the  counter  modulation  angle  deception  inter¬ 
ference  can  start  operating  is  that  the  radar's  useful  signals  are 
interrelated  on  the  wave  form  and  antenna  space  direction  angle. 
Previously,  when  we  deduced  the  radar's  integrated  resolution 
parameter,  it  was  assumed  that  the  signal  wave  form  and  antenna 
directional  angle  were  not  interrelated.  Because  of  this,  the 
obtained  radar  ECCM  capabilities  expression  did  not  include  this 
element  and  thus  awaits  further  research. 

V.  Summary  and  Calculation  Examples 

We  will  sum  up  in  table  2  the  radar's  ECCM  capabilities 
expression  proven  above. 
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Table  2 


Key: 


1*  Anti-passive  jamming  capabilities  (AC) 

2.  Anti-active  jamming  capabilities  (AJ) 

3.  Integrated  anti- jamming  capabilities  (AJC) 

4.  Fundamental  expression 

5.  Supplementary  factor 

6 .  Quick  changing  frequency  factor 

7 .  Side  lobe  factor 

8.  Moving  target  quality  factor 

9.  Constant  false  alarm  factor 

10.  Wide-limited-narrow  circuit  factor 

11.  Antenna  polarization  factor 

12.  Repetition  frequency  chatter  factor 


13. 

Whole  expression 

14. 

See 

formula 

(37) 

15. 

See 

formula 

(38) 

16. 

See 

formula 

(23) 

17. 

See 

formula 

(34) 

18. 

See 

formula 

(39) 

19. 

See 

formula 

(40) 

20. 

Or 

21. 

Or 

22. 

Or 

It  can  be  seen  that  the  radar's  ECCM  capabilities  expression 
has  two  parts.  The  first  is  the  fundamental  part.  This  part 
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measures  a  radar's  fundamental  or  potential  ECCM  capabilities 
determined  by  its  major  technical  parameters.  Its  significance 
is  similar  to  other  radar  fundamental  measurement  formulas  in 
wide  use  at  the  present.  For  example,  the  fundamental  measure¬ 
ment  formula  of  search  radar  force  is  the  power  and  antenna 
aperture  product  PA;  further,  the  fundamental  measurement  form¬ 
ula  for  the  accurate  tracking  radar 1 s  goniometric  angle  pre¬ 
cision  is  the  power,  antenna  aperture  and  (antenna  gain) [2]’ 
product  PAG [ 2 ]  [7], 


The  second  is  the  supplementary  factors  part.  This  part  shows 

the  radar's  ECCM  measures  and  their  quality  level.  Within  this 

part  those  factors  belonging  to  the  antenna  aspect  can  be  jointly 

called  antenna  quality  factor  FA  (only  the  part  that  acts  on  the 

passive  jamming  is  FA(J  .  Those  belonging  to  the  signal  aspect 

are  jointly  called  signal  quality  factor  F  (those  that  act  on 

s 

passive  jamming  are  Fs c  and  those  that  act  on  active  jamming  are 


By  combining  the  fundamental  part  and  supplementary  factors 
part  we  can  obtain:  anti-passive  jamming  capabilities  expression 
AC=T  B.GF..F  ;  anti-active  jamming  capabilities  expression 

AJ=PTqB^GFaFsj.  These  two  formulas  form  an  approximately  direct 

ratio  with  the  radar's  corresponding  output  signal  interference 
power  ratio.  In  the  formulas,  each  factor  cannot  compensate  for 
the  other. 
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Finally  we  have  the  combined  anti-jamming  capability  expression 
AJC  =  PTgB^GF^Fg -  This  formula  can  broadly  and  generally  indicate 

the  overall  capability  of  radar  to  resist  various  types  of  jamming 
but  is  is  not  directly  proportional  to  the  radar  signal  and  overall 
jamming  power  ratio.  In  the  formula,  the  various  factors  cannot 
compensate  for  each  other. 

Now,  we  will  discuss  the  measuring  method  for  a  multichan- 
neled  radar's  ECCM  capabilities. 

Many  radar  systems  consist  of  multiple  channels  and  they  are 
divided  into  several  different  categories: 

(1)  Each  channel  has  a  different  carrier  frequency,  yet  the 
space  beams  are  the  same  or  overlapping;  for  example,  frequency 
diversity  system  radar. 
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(2)  Each  channel's  carrier  frequency  is  identical,  yet 
the  space  beam  positions  are  different  or  only  have  partial 
overlapping.  For  example,  pulse  system  radar  and  some  three 
coordinate  radar. 

(3)  Each  channel's  carrier  frequency  is  different  and  the 
space  beam  positions  are  also  different  or  only  have  partial 
overlapping;  for  example,  some  three  coordinate  radar. 

For  each  channel's  ECCM  capability  measurement  formula,  as 
was  already  shown  above,  we  calculate  the  whole  radar  system's 
ECCM  capability  and  consider  them  separately  according  to  three 
categories  of  conditions. 

For  the  frequency  diversity  radar  which  belong  to  category 
(1) ,  one  channel  receives  interference  while  the  other  channels 
still  operate  normally.  If  all  the  channels  receive  interfer¬ 
ence  then  the  whole  system  is  inhibited.  Because  of  this,  the 
system's  anti-active  jamming  capacity  is: 

n 

wo,- 2  wo#  (40 

i  -  1 


In  the  formula,  n  is  the  channel's  total.  If  each  channel's 
carrier  frequency  differs  greatly  and  do  not  belong  to  the  same 
frequency  band  then  the  systems  anti-passive  jamming  capability 
is : 

mo,- 2  wo#  <<2> 

i»l 

Yet,  in  all  situations,  each  channel  carrier  frequency  still  be¬ 
longed  to  the  same  frequency  band.  At  this  time: 
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(AOt-lAOi, 


(43) 


In  the  formula,  (AC) .  represents  the  largest  one  in  each 
channel  (AC).  Using  the  system's  integrated  ECCM  capability 
calculation  method,  the  above  analysis  should  be  the  same  as 
the  anti-active  jamming  capability.  Therefore  it  is  defined  as: 


(AJC)s=  £  (MC), 

i«  1 


(44) 


For  the  category  (2)  multiple  beam  radar  system,  when  one 
space  beam  receives  interference  the  other  independent  beams 
still  operate.  Because  of  this,  the  system’s  anti-active  jam¬ 
ming  capability  is: 

II 

W,-i,J[(ij/)r  (45) 

■  » l 


In  the  formula,  kR  represents  the  overlapping  coefficient  of  the 
partial  overlapping  condition  in  the  beam.  If  each  beam  has  an 
independent  1^=1  and  if  there  is  a  partial  overlapping  of 
k  K  1  then  when  we  previously  defined  the  radar's  space  dir- 
ectional  angle  resolution  and  took  the  antenna  beam's  3 
decibel  angle  width  as  the  resolved  unit,  if  the  space  two  beams 
are  overlapping  below  3  decibels  this  can  be  viewed  as  the  fund¬ 
amental  independence,  taking  kR=l.  With  the  same  reasoning,  the 
system's  anti-passive  jamming  capability  is: 


(-JCWkS^Oi  (40) 

1*1 

The  system’s  integrated  ECCM  capability  is: 

n 

( AJC)s=kkV(AJC ),  0-) 

i  =  ! 

For  the  category  (3)  systems,  the  anti-active  jamming  capa 
bilities  are  the  same  as  formula  (41)  .  As  regards  the  anti¬ 
passive  jamming  capabilities,  we  observe  whether  or  not  each 
channel ' s  carrier  frequency  belongs  to  the  same  frequency  band 
and  use  formula  (42)  or  (46).  The  system's  integrated  ECCM 
capabilities  are  still  based  on  formula  (44) . 

Below  we  will  calculate  and  compare  four  types  of  radar 
ECCM  capabilities  according  to  each  of  the  above  mentioned  form 
ulas.  They  are:  two  new  models,  the  three  coordinate  radar 
AN/TPS-43E  (manufactured  by  Westinghouse  Company  of  the  U.S.) 
and  AR-3D  (manufactured  by  Bendix  Company  of  the  U.S.);  two 
older  search  radar  A  and  B  types  (see  table  3) . 
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Table  3. 


Key:  1.  Radar  system 

2 .  Remarks 

3 .  Type 

4.  Multiple  beam  6  channel  three  coordinate 
radar  AN/TPS-43E 

5.  Frequency  scanning  system  8  channel  three 
coordinate  radar  AR-3D 

6.  Multiple  beam  5  channel  two  coordinate 
radar  A  type 

7-  Frequency  diversity  system  4  channel  two 
coordinate  radar  B  type 

8 .  Radar  system  with  each  channel ' s  power 
and  antenna  gain  (see  table  4) 

9.  Total  emission  power  (kilowatts) 

10.  The  first  channel's  antenna  gain  G 

(decibels)  1 

11.  The  first  channel's  emission  power  P 

(kilowatts)  1 

12.  Emission  loss  (decibels) 

13.  Radar  signal  band  width  B  (megahertz) 

s 

14.  Signal's  continuous  time  T^  (seconds) 

15.  T  B  G.  (decibels) 

O  S  x 

16.  P.T  B  G.  (decibel  watts) 

l  o  s  l 

17.  Supplementary  factors:  k  (decibels) 

18.  (decibels)  J 

19.  (decibels) 

20.  (decibels) 

21.  (decibels) 

22.  (decibels) 

23.  (decibels) 

24.  (decibels) 

25.  (decibels) 

26.  Anti-passive  jamming  supplementary  factors 
total  (decibels) 

27.  Anti-passive  jamming  supplementary  factors 
total  (decibels) 

28.  Supplementary  factors  total  (decibels) 

29.  First  channel's  anti-passive  jamming  capa¬ 
bility  (AC)  (decibels) 

30.  First  channel's  anti-active  jamming  capa¬ 
bility  (AJ).,  (decibel  watts) 

31.  First  channels  integrated  ECCM  capability 
(AJC) ,  (decibel  watts) 

32.  System's  anti-passive  jamming  capability 

(AC)  ,  (decibels) 
s 
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33.  System's  anti-active  jamming  capability 

(AJ)  ,  (decibel  watts) 
s 

34.  System's  integrated  ECCM  capability 

(AJC)  ,  (decibel  watts) 
s 

35.  The  first  channels  of  the  first  three  types 
of  systems  are  the  beam's  lowest  elevated 
angle  channel 

36.  The  first  two  types  are  estimated  values 

37.  Estimated  value 

38.  Emission  loss  already  deducted 

39.  The  latter  three  types  of  radar  systems 
do  not  have  quick  changing  frequencies 

40.  All  are  non-polarized  adjustable  measures 

41.  B  type  without  MTI 

42.  All  without  Dicke-Fix 

43.  The  first  two  types  are  estimated  values  of 
the  A  and  B  types  which  are  without  CFAR 

44.  The  first  two  types  are  estimated  values 
of  the  A  and  B  types  which  are  without 
repetition  frequency  chatter  measures 
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Table  4  The  General  Characteristics  of  Each  Channel  for  Four 
Types  of  Radar 


Key:  1.  Name  of  radar 

2 .  A  type 

3 .  B  type 
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4.  The  operational  distance  for  each 

channel  is  R  ,  the  antenna  gain 
max 

is  G,  the  reflection  power  is  P  and 
there  is  a  mutual  proportion  rela¬ 
tionship 

5 .  Channel  number 

6 .  Other 

7.  Other 

8.  With  MTI 

9.  Without 

10.  Without 

11.  Without 

12.  Without 

13.  Without 

14.  Without 

15.  Without 

16.  Other 

17.  With  MTI 

18.  Without 

19.  Without 

20 .  Without 

21.  Other 

It  can  be  seen  from  the  calculation  results  that  because  the 
AR-3D  signal  band  width  is  large  and  its  antenna  gain  is  high, 
its  anti-passive  jamming  capability  is  about  6  decibels  higher 
than  that  of  the  TPS-43E.  Yet,  because  the  AR-3D  does  not  have 
quick  changing  frequency  capability,  its  anti-active  jamming 
capability  is  actually  about  4  decibels  lower  than  the  latter. 

It  is  worth  noting  that  the  AR-3D  systems  anti-passive  jamming 
capability  only  increases  0.1  decibels  as  compared  to  the  first 
channel.  Besides  the  first  channel,  this  is  due  to  the  other 
channels  not  having  MTI  circuits. 


As  regards  the  two  older  types  of  search  radar,  the  numer¬ 
ical  values  calculated  from  their  fundamental  ECCM  abilities 
expression  are  not  low  and  are  actually  close  to  those  of  the 
first  two  types  of  three  coordinate  radar.  This  shows  that  their 
ECCM  potential  capabilities  are  strong.  Yet,  because  they  bas¬ 
ically  do  not  add  any  ECCM  technical  measures  or  only  have  very 
few  measures,  the  supplementary  factors  are  all  relatively 
high  negative  values  and  the  final  calculated  ECCM  capability 
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expressed  numerical  values  are  much  lower  than  those  of  the 
three  coordinate  radar.  Their  performance  in  an  interference 
environment  drops  drastically  and  their  difference  from  modern¬ 
ized  radar  is  very  obvious.  On  the  other  hand,  in  view  of  their 
possessing  very  strong  ECCM  potential  capabilities,  if  we  carry 
out  technological  innovations  and  add  necessary  ECCM  technical 
measures  then  we  can  remake  it  so  that  it  possesses  ECCM 
capabilities  close  to  those  of  modernized  radar.  We  can  adopt 
formula  calculations  to  obtain  the  increased  numerical  value 
index  of  their  ECCM  capabilities. 

Now  we  can  sum  up  the  major  uses  of  the  radar’s  ECCM  capa¬ 
bilities  measurement  formula  which  we  obtained  as  follows: 

(1)  We  can  use  the  radar's  own  technical  parameters  to 
measure  the  numerical  value  indices  of  the  radar's  anti-passive 
jamming  capabilities,  anti-active  jamming  capabilities  and  in¬ 
tegrated  ECCM  capabilities.  This  is  convenient  for  making  com¬ 
parisons  between  various  types  of  radar  or  between  various  types 
of  design  plans  for  certain  radar. 

(2)  Because  the  anti-passive  jamming  capability  expression 
and  anti-active  jamming  capability  expression  are  approximately 
in  direct  ratio  to  the  output  signal's  interference  power  ratio 
under  corresponding  interference  conditions,  we  can  take  their 
calculated  numerical  values  and  directly  use  them  to  compare 
and  calculate  the  levels  of  surveyed  location  performance 
decreases  for  the  various  plans  of  various  types  of  radar  or 
one  type  of  radar  in  an  interference  environment. 

(3)  The  radar's  ECCM  capabilities  expression  includes  the 
fundamental  and  supplementary  factors  parts.  Because  of  this, 
we  can  separately  calculate  the  radar's  potential  ECCM  capabil¬ 
ities  and  the  effects  of  the  radar's  ECCM  technical  measures. 

For  certain  radar  we  can  use  this  analysis  in  design  or  in 
technological  innovation  to  increase  the  investment-result 
ratio  of  certain  complex  ECCM  measures. 
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Finally,  it  should  be  pointed  out  that  when  the  same  radar 
is  in  different  operational  conditions  (for  example,  search  or 
tracking  conditions) ,  if  its  technical  parameters  or  channel 
numbers  are  different  then  the  ECCM  capabilities  numerical 
values  calculated  according  to  the  formulas  will  be  different. 


References 

1.  S.L.  Johnson,  ECCM  Improvement  Factors (EIF) ,  International 
Countermeasures  Conference,  1976. 

2.  Li  Nengjing:  Measuring  Method  for  Radar's  ECCM  Capabilities, 
Acta  Aeronautica  et  Astronautics  Sinica,  no.  1  (1979) , 

pp.  70-81. 

3.  H.  Urkowitz,  C.A.  Hauer  and  J.F.  Koval,  Generalized  Resolution 
in  Radar  Systems,  Proc.IRE.  50,  No.  10,  Oct.  1962, 

pp.  2093-2105. 

4.  A.W.  Lihaijieke  (?  roman izat ion ) :  The  Radar  Resolution  Theory, 
Scientific  Publishers,  1973. 

5.  Zhang  Zhizhong:  The  Selection  and  Processing  of  Radar  Signals, 
Defense  Industries  Publishers,  1979. 

6.  B.D.  Steinberg,  Signal  Enhancement  by  Linear  Filterinq  in 
Pulse  Radar,  Lectures  20-22  of  special  summer  course, 
University  of  Pennsylvania  (June,  1961) . 

7.  D.K.  Barton,  Real-World  Radar  Technology,  The  Record  of  the 
IEEE  1975  International  Radar  Conference,  pp.  1-22. 


181 


Abstract 


This  article  reviews  the  commonly  used  methods  for  expression  of  radar's 
ECCM  performance*  points  out  their  disadvantages  and  proposes  that  it  is  re¬ 
quired  to  establish  the  formulae  which  consist  only  of  the  radar' s  technical 
parameters  in  order  to  measure  the  radar's  ECCM  capabilities. 

By  means  of  logical  argumentation  and  ambiguity  function  theory  we  can 
derive  the  followings*  (l)  The  radar's  4-dimensional  generalized  resolution 
parameter  can  be  used  to  express  the  radar's  basis  anti-clutter  capability.  Thf 
simplified  expression  of  this  parameter  is  T,  (the  dwelling  time  at  a  radar's 
target,  or  the  integration  time  of  a  echo  signal  of  the  radar)  x  Ba  (the  radar 
signal  bandwidth)  X  G  (the  radar  antenna  gain)i  (  2  )  The  product  of  an  ave¬ 
rage  radar's  transmitting  power  and  its  generalized  resolution  parameter  PT0  x 
BsG  can  be  used  to  express  the  radar' s  anti-jamming  capability,  and  the  sa¬ 
me  formula  may  also  be  used  to  express  the  total  interference  rejection  capa¬ 
bility  in  gross. 

By  another  approach,  i.  e.  by  calculation  of  the  radar's  output  signal  to 
interference  power  ratio  under  various  interference  conditions,  the  given  for¬ 
mulae  have  further  been  proven. 

The  given  formulae  consisting  of  the  radar's  fundamental  parameters  are 
the  fundamental  parts  of  *he  measuring  formula  of  the  radar' s  ECCM  capuhi- 
lit'es.  They  represent  the  radar  basic  or  potential  ECCM  performance.  The 
supplementary  factors  for  the  ECCM  capability  formulae  are  composed  of  the 
technical  specifications  of  the  radar's  ECCM  devices.  The  general  rule  of  com¬ 
position  is  to  take  the  typical  quality  level  (or  the  middle  level)  specifica¬ 
tion  (in  db)  as  the  base  figure.  If  the  quality  level  of  some  devices  is  higher, 

a  positive  factor  will  be  got,  and  if  lowei - a  negative.  Cor  example,  the 

equation  for  MTI  supplementary  factor  Pu  is 

P*(db)  =  SCV  — 25(db) 

where  25  db  is  the  middle  level  specification  for  SCV  of  modern  MTI  radar. 

As  many  radar  systems  consist  of  multiple  channels,  in  this  article  we 
propose  a  set  of  equations  for  the  multichannel  radar  system's  ECCM  capa¬ 
bilities  in  terms  of  the  ECCM  capabilities  of  the  individual  channel. 

In  the  last  part  of  this  article,  4  radar  systems  are  taken  as  examples, 
their  ECCM  capabilities  are  calculated  by  means  of  the  given  formulae.  Appli¬ 
cation  of  obtained  results  to  radar's  system  analysis  and  system  design  are 
discussed. 
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BRIEF  REPORT  ON  THE  FIRST  FLIGHT  TEST  RESEARCH  ACADEMIC 
EXCHANGE  MEETING 


The  Chinese  Aeronautics  and  Astronautics  Institute  con¬ 
vened  China's  first  flight  test  research  academic  exchange 
meeting  from  December  16-21,  1980  at  Yentai  in  Shandong 
Province.  The  conference  was  prepared  and  managed  by  the  Test 
Flight  Institute.  Attending  the  conference  were  28  representa¬ 
tives  from  related  units  and  77  authors  of  papers  from  China. 
Among  these,  besides  the  engineers  and  technicians,  there 
were  also  pilots  with  long  term  test  flight  experience. 

At  the  conference,  74  papers  were  exchanged.  The  contents 
were  concerned  with  the  flight  test  theories,  methods,  test 
techniques,  flying  techniques  and  testing  techniques  of  air¬ 
craft,  dynamic  devices,  electronics,  special  designed  air  con¬ 
ditioning  and  aviation  armaments.  The  conference  was  divided 
into  three  stages.  The  first  stage  was  a  general  meeting  ex¬ 
change  and  besides  exchanging  seven  papers,  two  special  reports 
were  given.  The  second  stage  was  divided  into  three  sets  of 
exchanges:  aircraft,  dynamic  devices,  and  electronic,  special 
designs  and  test  instruments.  Altogether,  44  papers  were  read 
(another  23  written  papers  were  exchanged) .  The  third  stage 
consisted  of  specialized  discussions  of  problems  dealing  with 
the  advance  of  flight  test  research,  promoting  academic 
exchange  and  establishing  related  academic  organizations.  From 
beginning  to  end  the  conference  was  very  conscientious,  enthus- 
iatic  and  pervaded  with  a  friendly  atmosphere. 
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The  conference  had  abundant  content,  wide  specialized 
fields  and  strong  real  practice  which  caused  everyone  to 
broaden  their  outlook,  increase  their  knowledge  and  attain  to 
great  gains.  Everyone  considered  that  this  indicated  China's 
flight  test  research  enterprises  were  going  from  the  stage  of 
having  studied  and  comprehended  foreign  test  flight  theory  and 
methods  and  were  beginning  to  enter  a  stage  in  which  China 
was  initiating  its  own  independent  new  creative  research.  The 
representatives  each  requested  that  this  type  of  academic 
meeting  be  convened  in  the  future. 

During  the  conference,  they  discussed  existing  problems  in 
China  's  present  test  flight  research,  proposed  establishing 
corresponding  academic  organizations  and  chose  a  name  list  for 
the  preparatory  group  of  the  academic  organizations.  The  con¬ 
vening  of  the  second  academic  exchange  meeting  during  the  third 
quarter  of  next  year  by  the  Flight  Mechanics  Alliance  was  also 
decided  at  the  conference. 
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AIMING  COMPUTATION  FOR  FIGHTER  WEAPON  AIMING  SYSTEM 
Zhang  Sen 

(Optic  Machinery  Research  Institute) 

Abstract 

This  paper  presents  several  derivations  of  the  aiming 
computing  equations  for  the  "Non-Director"  (also  called  "Dis¬ 
turbed  Recticle")  system  in  the  aiming  computer  systems  of  the 
air-to-air  gunsight  and  head-up  display  weapon.  This  paper 
deals  with  four  situations:  the  fighter  attacks  a  non-maneuver- 
able  target  which  is  considered  as  moving  in  a  straight  line 
at  a  constant  speed;  the  fighter  attacks  a  maneuvering  target 
and  the  corrections  for  the  maneuvering  of  the  target  are  taken 
into  account  in  the  aiming  computation;  the  fighter  attacks  a 
maneuvering  target  and  the  corrections  for  the  maneuvering  of 
the  target  and  the  effect  of  the  fighter's  own  roll  rate  are 
also  taken  into  account  in  the  aiming  computation;  during  the 
attack  of  the  fighter  on  a  maneuvering  target  the  historical 
tracer  line  (hot  line)  and  a  version  of  the  damped  tracer  line 
based  on  multiple  Lcos  computations  are  adapted  in  th«.  com¬ 
putation  for  air-to-air  gun  snap-shot. 

Symbols 

M  target  point 
My  target  leading  point 

%  bullet's  oblique  firing  length 

n  ballistic  lowering  quantity 

5  target  distance 


D,D 


V01 

V° 

V01 


*M 

7 

cp 

U) 

M 


°-XlylZl 


a,  0 
0 

r 

Vim. 

t*.v(~v  } 


Mxx  ,  WMy^  A>Mzx 

^Ixr  Vr  6)1  z1 

Tj 

Ty 


distance  change  rate  and  distance  change 
acceleration 

commencing  speed  of  muzzle 
fighter's  speed 

bullet's  composite  velocity  (vgi=V0  +  Vl^ 

unit  vector  of  Vq^ 

target  velocity 

bullet's  mean  velocity 

target's  relative  angle  velocity 

fighter's  rotational  angle  velocity 

is  the  fighter's  body  coordinate  system 
xlylzl  are  separately  the  unit 

vectors  on  three  axes 

fighter’s  attack  angle  and  slip  angle 

fighter's  pitching  angle 

fighter's  slope  angle 

pitching  total  correction  angle  (or 
pitching  total  leading  angle) 

azimuth  total  correction  angle  (or 
azimuth  total  leading  angle) ;  when  u  and  v 
are  in  air-to-air  gun  snap-shot,  this  also 
indicates  the  fighter's  pitching  and 
azimuth  angles  of  the  tracer  point  on  the 
tracer  line. 

indicate  the  projection  of  the  target's 
relative  angle  velocity  on  the  three  axes 
of  the  fighter  body's  coordinate  system 

indicate  the  projection  of  the  fighter's 
rotational  angle  velocity  on  the  three 
axes  of  the  fighter  body's  coordinate 
system 

*  ~ 

computation  time  (Tj=aTj) 

the  flight  time  required  from  after  the 
bullet  is  shot  until  it  hits  the  target 

the  rate  of  change  for  the  pitching  and 
azimuth  total  correction  angles 
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^pitching'  ^azimuth 

n 

alift 

a tangential ' 
acentripal 
S 

^1 

i 

ab,  Klt  K2,  K3,  F,  a 


the  fighter's  pitching  and  azimuth 
aiming  angle  components 

fighter's  normal  overload 

fighter's  lift  acceleration 

the  fighter's  tangential  and  centripal 
accelerations 

Laplace  operator 
time  constant 

gravitational  acceleration 

the  constants  related  to  the  ballistic 
and  fighter  flight  parameters 


All  of  the  vectors  have  the  following  relationship  with  the 
fighter  body's  coordinate  system: 


d=*D Ceos n cos v  sinH  —  cost* sin  v  j 


D=>b  CcoaUcosv  sinU  —  cos  H  sin  v] 


,  (cos  a  cos  0  sina  —  cosasinPJ 


5#»a#C  0  1 

(1)  (2) 


0  3 


H  C— sin  0  — cos  0  cosY  cos  9  sin  Y  3 


Yi 

-  *1  J 

.  K 
*»’ 

*1  j 
*1 

91 

91 

*»  J 


Key:  1.  Lift 

2.  Lift 
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Fig.  1.  The  Relationship  of  All  of  the  Vectors  and  the  Fighter's 
Fuselage  Coordinate  System 

In  the  fighter's  fuselage  coordinate  system  o-x^z^ 

0  is  the  point  of  origin  in  the  fighter's  gravi¬ 
tational  center?  the  ox^  points  straight  toward 
the  nose  along  the  aircraft's  longitudinal  axis 
direction;  the  5YT  points  straight  toward  the 
cockpit  along  theiaircraf t ' s  mechanical  axis 
direction;  and  ojT^  ,  and  oy^  form  a  right  angle 

right  hand  coordinate  system  relationship. 

Key;  1.  Lift 


The  total  design  of  the  fighter  weapon  aiming  system  has  a 
multifaceted  content  and  the  determining  of  the  aiming  computer 
equations  is  one  of  its  important  functions.  Below  we  will  dis¬ 
cuss  this  problem  under  four  types  of  different  conditions. 

1.  The  aiming  computation  of  a  non-maneuverable  (or  rela¬ 
tively  lacking  maneuverability)  air  target. 
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Fig.  2  Aiming  Vector  Diagram 

Fig.  2  is  the  aiming  vector  diagram  under  this  type  of 
condition  and  from  this  we  can  arrange  a  deviation  vector  aiming 
equation : 

N  =  D+VuT,-lyl{-  ij 

Because 


therefore 


V«=Vi  +  -^*=  V,  +  D+(a,,x5) 


► 


N  =  D+  VtT,+  or,+  (%XD)r,-  4Vo,-ij 

=(i  +  g- T,)  d+(v»xd)t,+v,t,- 
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Because 


4 


therefore  we  obtain 


1  +~d~t*)  d+(uuxd) 


Given  that  the  gun  axis  and  fuselage  axis  coincide,  we  use 
the  relationship  of  all  of  the  vectors  and  the  fuselage  coordinate 
system  and  project  it  toward  the  fighter's  fuselage  coodinate 
system.  Thetefore,  a.P.,u.,  v,  *1  and  sin  0  are  all  small  and 
relatively  small  quantities  and  because  of  this  we  can  consider 
sin  a**a,  sin  P  **  0.  sin  uzsu,  sin  vjsv,  cos  ^czcos  ^xcoswsscos  v 
css.  1 .  If  we  overlook  the  quadratic  or  above  quadratic  small 
quantity  then  we  can  obtain  the  computing  equations  for  (N=0) 
pitching  channel  total  correction  angle  Ptet*)  and  azimuth 
channel  total  correction  angle  ( Tfe.v  )  under  aiming  conditions : 


^»i»= 


■00 


k.,7'/  + 


'l  cos  0  cos  Y 


„  ~  *1  cos  8  sin  Y 

,<W'-(f.;-p)7'7 


In  the  equations,  Tj=  ^ — ^tt-  is  the  smoothing  of  the  aiming  line 

vcp  1 

during  track  aiming  and  this  causes  the  return  of  the  peonle  and 
aircraft  to  be  steady.  Damping  coefficient  x  is  introduced  into 
the  aiming  computing  equation;  to  obtain  the  target's  relative 
angle  velocity,  we  introduce  the  relationship  of  the  angular 
velocity  and  the  target's  relative  angle  velocity  when  the  fighter 
is  track  aiming:  a>,  =  +  "ts.  r  to  decrease  the  dynamic  error  during 

track  aiming  we  take  T^=3Tj  (5  is  smaller  than  the  coeffient  of 
1) .  In  this  way  we  obtain  the  following  equations: 
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,  v .  ,  1  cos  ft  cos  Y 

:  (tv-r.rr, 

/  ,  v  -  ,  V,  ’l  cos  0  sin  Y_ 

+*,-C®i»,-(  1  +  x  KtO  1 1~  (F^-VjT, 


For  the  sake  of  simplification  and  convenience,  if  we  can 
consider  that  during  an  attack  the  fighter  makes  equal  overload 
circles,  then  the  second  item  (aiming  angle  item)  on  the  right 
end  of  the  equation  can  obtain  the  following  form: 


Azimuth  aiming  angle 
Key:  1.  Azimuth 


_  *1  cos  9  sin  Y  _  *1  cos  0  J'rlcolcosY 

UT(<  0  cos  0 

_  V*»  V„Wl 

’ii\,-v,)ngnr~D,of£ 


Pitch  aiming  angle 
Key:  1.  Pitch 


_  n  cos  0  cosy  _ 

(i)  (v„-vx)T,  D,ar 


cos  0  cos  Y 


To  correct  the  azimuth  aiming  angle  then  we  use  the  method 
of  multiplying  Tj  times  a  coefficient  K2  smaller  than  1  to 
approximate  the  correction  and  obtain: 


C®i»it~(  1  +  *  H*07  *Kt 


The  total  correction  angle  of  the  pitching  channel  also  changes 
because  of  this: 


G) 


l»I 


V  n  T* 

COSY 


It  is  simplified  into: 


V  n  r*  f 

'C®i,,-(  1  +  *  \ 


cos  9  cos  Y  +- 
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Because 


(aT*  +  6 ) 


■cos  9  cos Y  +  - 

Key:  1.  and  take 


we  then  obtain:  1  +  «  );*»)T*Kl  +  (aH+  6  )  " 

Naturally,  in  accordance  with  the  difference  of  the  selec¬ 
tion  and  simplification  methods  for  the  aiming  system  plan,  we 
can  also  obtain  other  forms  of  the  aiming  computing  equation. 
For  example: 


—  (  1  +  *  +  (aT*  +•  b  ) 

'|it>  ~  C®i»|  “  (  1  +  *  )  \(ix»3 1  *Ki 


Whether  it  is  the  aiming  computing  accuracy  or  decreasing 
the  system's  complex  equations  both  are  satisfactory. 

During  the  previous  derivation  and  simplification  process 
for  the  aiming  computing  equations,  we  overlooked  the  influence 
of  the  trailing  angle  item  when  shooting;  further,  because  the 
attack  is  against  a  non-maneuverable  (moving  in  a  straight  line 
at  a  constant  speed)  target,  the  roll  rate  (the  angular  rate 
projected  on  the  fuselage  coordinate  system  x^axis)  generally 
does  not  exceed  3°/second  and  the  error  caused  by  it  is  usually 
within  2-3  milliradian.  Thus,  we  do  not  have  to  consider  it. 

2.  The  approximate  aiming  computing  correction  for  a 
maneuvering  target 

This  type  of  aiming  computing  correction  is  for  improving 
the  original  sight  (non-head-up  display)  weapon  aiming  system 
for  use.  Its  aiming  vector  diagram  is  shown  in  figure  3. 
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Fig-  3.  Aiming  Vector  Diagram 

The  arrangement  of  its  deviation  vector  aiming  equation  is: 

N*DtFkT,+  f  o 

When  the  bullet  flight  time  is  relatively  short  (close  to 
the  range  of  fire) ,  we  take  VM  as  a  constant  (equal  to  VM  when 
opening  fire) . 

Further,  Vi  +  5  +  («**  d) 

p„-pl+2)  +  2miix5+«mX  («mXd)-+  Chxd) 


Overlooking  the  angle  and  second  derivative  item  and  second  small 
quantity  of  the  distance,  we  then  obtain: 

N-d+v,T,+  dT,+  (u,,x  5)  Tr+jvin+(»«XD)71--5v;i-ii 


192 


When  the  fighter  is  in  attack  flight,  V,  =a,_  .  +  a  .  , 

^  ^  1  tangent  vector' 

yet  ^tangent  has  a  very  sina11  influence  on  the  aiming  computation. 
Therefore,  we  can  take  V ^  avector;  take  n^^gT^  (from  this,  the 
caused  error  is  generally  within  1*) ;  when  the  fighter  is  in  equal 
overload  circles,  there  is  the  following  relation  in  the  accelera¬ 
tion:  ilift=itangent  +  (-5) •  Thus  we  can  obtain: 

,  (1)  . 

N  -  D  +  V  ,7V +  dT, +  (au  x  d)  T,  +  ya»n  +  (S*  x  5) 71 - 1 vlt 
Key:  1.  Lift 

Projected  on  the  fighter's  fuselage  coordinate  system,  under  aim¬ 
ing  conditions,  after  N=0  simplification,  we  obtain: 

Key:  1.  Lift 
In  the  formula,  a»=*ng 

9 _ 

K'~  2{V.,-V,) 

t  +  x  )*tJ7y 
Key  1:  Lift 

This  aiming  computing  equation  took  into  account  the 
correction  of  the  maneuvering  target  but  did  not  take  into  account 
the  influence  of  the  roll  rate.  This  is  because  most  sight  com¬ 
puters  are  electromechanical  or  electronically  simulated  wherein 
it  is  difficult  to  realize  the  correction  of  the  roll  rate 
influence.  Yet,  it  must  be  pointed  out  that  when  there  is  a 
relatively  large  roll  rate  (over  10°/second)  the  caused  error 
cannot  be  overlooked.  For  example,  when  a  fighter  carries  out 
aiming  attack  on  a  target  with  equal  overload  circles,  the 
application  of  this  formula  group  for  aiming  computing  will  be 
very  accurate.  Furthermore,  this  set  of  equations  can  also 
realize  the  computing  of  the  reticle  snap-shots  which  does  not 
require  a  central  spot  of  the  reticle.  It  does,  however,  require 
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that  before  the  reticle's  central  spot  goes  on  the  target  we 
use  the  method  of  advancing  the  flight  time  firing  of  a  shell 
and  only  then  will  the  target  be  brought  down. 

3.  Weapon  aiming  computing  taking  into  account  the  correc¬ 
tions  for  the  maneuvering  of  a  target  and  the  influence  of 
the  fighter's  roll  rate. 

Its  relationship  within  the  aiming  vector  can  be  seen  in 
fig.  3.  The  deviation  vector  aiming  equation  is: 

N=D  +  V-ry+  |*  *  Vudldi-lvls-v 

With  the  same  principle  as  was  mentioned  previously  in  the  second 

part,  V..  is  viewed  as  a  constant  and  when  the  relation  of  V,,  and 
M  M 

VM  is  substituted  into  the  equation  and  we  overlook  the  angle 

and  second  derivation  item  of  the  distance,  then  we  obtain: 

N-Zi  +  Vi  7’,  +  d7’,+  (w«xo)  T,  +  YViTl+(auxp)  TJ+ 
x  (a„x  D)Tl-lvh-n 

With  the  same  rationale  we  took 

Vi*''®*,  ’i*5  =  + 

(1)  (2)  (3) 

Key :  1 .  Tangent 

2.  Lift 

3 .  Tangent 

and  then  we  obtained 


N  -  D  +  V  J,+  dT,  +  (a*  x  o)T,  +  -~a*  7*  +  (a„  X  5)  T\ 

+  j  “M  (««x6)n-4vli 

Key:  1.  Lift 

In  the  formula 

a*  x  (»*  x  5)  -5<*>J 

When  this  formula  is  substituted  into  the  N  relation,  we  obtain: 
N“D+V17t,+57\,+  (S„xd)7\  +ya*  T;+(««x5)TJ 

£k»„in— l-DmiTj-lyJ, 

Key:  1.  Lift 

When  each  vector  is  projected  toward  the  fuselage  coordinate 
system  and  simplified,  we  can  obtain  the  aiming  computing 
equations  under  aiming  conditions: 

*«■—  v  i  +  *  i-cD,.,!!?, 

Key:  1.  Lift 

In  the  equations,  alift=ng. 

This  takes  into  account  the  correction  for  the  maneuvering 
of  the  target  as  well  as  the  fuselage ' s  roll  rate  influence  on 
the  corrected  aiming  computing  equation.  It  is  generally  applied 
for  use  in  the  numerical  computer's  weapon  aiming  system. 
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4.  Its  application  in  the  head-up  display  aiming  system  can 
carry  out  snap-shot  tracer  line  aiming  computing  for  air-supported 
maneuvering  targets. 

The  typical  display  of  this  type  of  aiming  attack  method 
developed  during  the  1970's  is  shown  in  fig.  4. 


Fig.  4.  Typical  Display  Diagram  of  Tracer  Line 

Key:  1.  Gun  cross  line 

2 .  0.5  seconds 

3.  Tracer  line 

4.  Reticle 

5.  Distance  octagon 

6.  1.5  seconds 

7 .  1  second 

8 .  Central  spot 

9.  Target  aircraft 


In  principle,  the 'snap-shot  tracer  line  can  be  divided  into 
two  types:  one  type  is  the  past  historical  tracer  line  (also 
called  the  hot  line)  and  the  other  type  is  the  damped  calculated 
tracer  line.  The  so-called  historical  tracer  line  indicates  the 
connecting  line  of  the  real  time  calculation  and  each  different 


time  a  hypothetical  bullet  is  fired  within  a  past  period  of  time 
opposite  a  fighter's  bullet  tracing  point.  Its  principle-vector 
relation  diagram  is  shown  in  fig.  5. 


Fig.  5.  Principle-Vector  Relation  Diagram 

Key:  1.  Firing  point 

2.  Bullet  point 

3.  Fighter's  realized  position 

4.  Fighter's  locus 


In  the  figure,  o  is  the  reference  point,  xoy  is  the  se¬ 
lected  coordinate  system;  A(t)  is  the  target’s  maneuvering 
acceleration. 


The  realized  bullet  point  position  is 


b  =  a«+4vJ  ,  +  »i 


The  realized  fighter  position  is 
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A-A.  +  J  #  J  #A(  i  )didi+vlT, 

Because  R=B-A,  R  is  the  vector  of  the  realized  bullet  position 
opposite  the  fighter's  realized  position.  Therefore,  we  obtain 

£=4  vl,+n  -f  o'f  o  A<t)did:+rtT, 

A(t)  is  the  fighter's  motion  acceleration. 

when  there  are  projection  calculations  of  all  of  the 
vectors  in  this  vector  equation  within  the  realized  fuselage 
coordinate  system  we  can  then  obtain  the  aiming  computing 
equation  of  the  past  historical  tracer  line.  Using  the  equation, 
we  can  calculate  the  bullet 1 s  tracer  point  each  time  a  hypo¬ 
thetical  bullet  is  fired  opposite  the  fighter’s  realized  position 
The  connecting  line  of  these  points  is  displayed  on  the  head-up 
display  which  is  the  past  historical  tracer  line. 

The  damped  calculated  tracer  line  is  sometimes  also  called 
the  damped  tracer  line.  Because  of  the  differences  of  the  com¬ 
puting  method  and  damped  lead-in  method,  it  has  many  different 
types.  The  most  commonly  used  among  them  is  the  type  of  tracer 
line  obtained  after  repeated  calculations  using  the  total 
correction  angle  aiming  computing  equations  derived  in  part 
three  above.  The  damping  in  the  equation  can  set  up  a  composi¬ 
tion  that  is  changeable.  The  total  correction  angle  aiming 
computing  equation  mentioned  previously  can  also  be  rewritten 
in  the  following  form: 


i+ls['""+  2  (F^rr  +  -f“"v]7'’ 

T+^s  [(0'> 

In  the  equations,  u  and  v  separately  indicate  each  different 
Ty  time  tracer  point  on  the  tracer  line  opposite  the  fighter's 
pitching  angle  and  azimuth  angle? 

1  is  the  damping  time  constant  which  can 
change  within  0.5  T,- 1.20  T,  • 

S  is  the  Laplace  operator. 

The  above  can  be  used  for  the  tracer  line  of  air-to-air 
gun  snap-shot  and  has  often  been  applied  to  antiaircraft  guns 
(such  as  the  American  six  barrel  gun) .  Furthermore,  if  this 
type  of  attack  method  is  used  for  most  single  barrel  aerial  guns 
with  relatively  low  firing  speeds,  the  shooting  down  probability 
of  air-supported  targets  cannot  be  very  high.  Naturally,  if 
several  more  simultaneous  volleys  are  installed  on  the  aircraft, 
in  principle,  this  is  also  a  type  of  remedy  method  yet  this  is 
further  limited  by  other  conditions. 

This  paper  only  analyzed  and  discussed  several  specific 
problems  regarding  the  fighter  weapon  aiming  system.  It  did  not, 
however,  mention  the  director  type  weapon  aiming  system  or  the 
MRGS  system  which  is  now  being  developed  abroad. 
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Abstract 

This  paper  presents  several  derivations  of  the  aiming  computing  equations 
for  the  “Non-Director”  or  the  “Disturbed  Reticle”  system  in  the  aiming  compu¬ 
ter  systems  of  the  air-to-air  gunsight  and  head-up  display  weapon.  The  paper 
deals  with  four  casest  the  fighter  attacks  the  non-maneuverable  target  which 
is  considered  as  moving  in  a  straight  line  at  a  constant  speedi  the  fighter  attacks 
a  maneuvering  target  and  the  corrections  for  the  maneuvering  of  the  target 
are  taken  into  account  in  the  aiming  computation!  the  fighter  attacks  a  ma¬ 
neuvering  target  and  the  corrections  for  the  maneuvering  of  the  target  and  the 
effect  of  fighter's  own  roll  rate  are  taken  into  account  in  the  aiming  compu¬ 
tation!  during  the  attack  of  the  fighter  on  a  maneuvering  target  the  historical 
tracer  line  (hot  line)  and  a  version  of  the  damped  tracer  line  based  on  mul¬ 
tiple  Lcos  computations  are  adapted  in  the  computation  for  air-to-air  gun 
snap-shot. 
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THE  OPTIMUM  DESIGN  OF  NON-MOMENT  LAMINATED  COMPOSITE  PLATE 
-  According  to  static  failure  strength  condition  - 

Ma  Zukang 

(Northwestern  Polytechnical  University) 

Abstract 

When  the  upper  and/or  lower  panels  of  the  loading  box  of 
the  aircraft  wing  surface  is  made  of  the  fiber  reinforced  com¬ 
posite  laminates,  they  can  frequently  be  simplified  as  a  non¬ 
moment  plate.  This  paper  introduces  an  optimum  design  method 
of  the  laminated  plate  according  to  static  failure  strength 
condition.  The  mathematical  tool  used  in  the  measure  is  the 
Lagrangian  multiplier  method  and  the  static  failure  strength 
condition  is  adopted  as  Hill-Tsai  criteria  and  Norris  criteria. 

The  formulas  for  optimum  design  have  not  only  been  derived, 
but  also  reformed  to  be  convenient  for  application  in  the  computer. 
How  to  establish  the  ultimate  strength  of  luminates  is  also  dis¬ 
cussed  in  detail.  An  example,  illustrating  the  solution  procedure 
and  how  to  select  the  optimum  scheme  of  lamination  design,  is 
presented  in  this  paper.  Some  technical  problems  are  briefly 
discussed  in  the  last  part. 

At  the  stage  of  the  preliminary  structural  design,  this  pro¬ 
cedure  can  be  considered  as  an  engineering  method  of  lamination 
optimum  design  for  the  loading  panel  of  laminated  composite 
which  works  under  tension  or  compression  (assuming  that  the 
buckling  failure  would  not  occur) . 


I.  Preface 

The  total  bearing  action  of  the  upper  and/or  lower  panels 
of  the  loading  box  of  the  aircraft  wing  surface  is  mainly  the 
axial  load  and  shear  load  in  the  bearing  plate  plane.  Even  the 
aerodynamic  load  of  the  vertical  plate  plane  is  a  local  bearing 
problem  which  cannot  be  overlooked  in  the  preliminary  structural 
design  when  calculating  the  total  strength.  Aside  from  this,  to 
avoid  the  coupling  effects  of  the  laminated  composite  plate,  the 
lamination  should  be  designed  into  a  symmetrical  mirror  surface. 
Furthermore,  because  the  curve  of  the  wing's  stressed  panel  is 
very  small  and  is  close  to  a  flat  plate,  the  wing's  stressed 
panel  can  be  taken  as  a  non-moment  laminated  plate.  That  is, 
moments  without  internal  or  external  force.  Further,  the  lamin¬ 
ated  plate  is  in  a  plane  stress  condition.  Frequently,  a 
0°/i45°/90o  lamination  is  used  by  this  type  of  loading  laminated 
plate  (the  principle  of  this  method  is  still  appropriate  for 
other  lamination  angles) . 

Because  the  lamination  of  the  laminated  plate  mainly 
affects  the  layer's  normal  stress  and  shear  stress  but  has  a  very 
small  effect  on  the  bearing  capability,  the  main  content  of  this 
paper  concerning  the  optimum  design  of  non-moment  laminated 
plate  is  the  determination  of  the  number  of  layers  needed  by 
each  directional  layer  based  on  the  size  of  the  external  load 
so  that  the  laminated  plate  will  have  the  lightest  weight  under 
certain  constraint  conditions.  These  constraint  conditions  are 
under  the  action  of  the  design  load  and  in  accordance  with  the 
static  failure  strength  criteria,  the  laminated  plate's  safety 
tolerance  is  zero  (or  the  remaining  strength  coefficient  is  1) . 
Mathematically,  this  is  the  problem  of  finding  the  conditional 
extreme  value  which  can  be  solved  by  using  the  Lagrangian 
multiplier  method. 

II.  Constraint  Conditions 

We  can  select  suitable  failure  criteria  based  on  the  type 
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a 


gx-.-a-Mg, 


reinforced  fiber.  For  example,  we  can  use  the  Norris  criteria 
(generalized  divalent  interaction  law) 


or  the  Hill-Tsai  criteria 


In  the  formulas,  the  N^,N^  and  are  the  design  loads  (service 

loads)  of  the  unit  lengths  on  the  x,y  and  s  axes  defined  in 
figure  1.  The  unit  is  pounds/inch  or  kilograms/centimeter. 


N , 


Fig.  1.  The  Plane  Stress  Action  on  the  Laminated  Plate 

In  N.=fvt=f.  (t  +t  =t  ),  (k  is  the  service  stress  of  the 
K  K  K  x  y  s 

k  direction  and  k=x,y,s. 

D  , D  and  D  are  the  ultimate  loads  (permissable  loads)  of 
x  y  s 

the  unit  lengths  on  the  y  and  s  directions.  The  unit  is  also 


pounds/ inch  or  kilograms/centimeter. 

In  Dk=F^t=F^(tx+ty+tg) ,  F^  is  the  ultimate  stress  (ultim¬ 
ate  strength)  of  the  k  direction  and  k=x,y,s. 

III.  Determining  the  Ultimate  Strength  of  the  Laminated  Plate 

When  we  take  into  account  that  the  ultimate  strength  of  the 
laminated  plate  along  a  certain  axis  direction  is  obtained 
from  the  ultimate  strength  of  each  directional  layer  provided  in 
the  corresponding  direction  according  to  the  linear  weighted 
average  rule,  it  is 

7  -j~  +  t7  “^+'F7  T")F* 

”  (a„Z.  +  a,yM  4-  a.,N)  F. 


In  the  formula, 

F  is  the  ultimate  strength  of  the  lamin- 
x  ated  plate  along  the  x  axis  (the  unit 

2 

is  kilograms/centimeter  or  1000 
2 

pounds/inch  )  ,* 

F  F  F 

xx '  xy'  xs  are  the  ultimate  strengths  able  to  be 
provided  along  the  x  axis  by  each  dir¬ 
ectional  layer  under  plane  stress 
conditions ; 

F  is  the  longitudinal  tensile  ultimate 
a  stress  of  the  single  direction  compos¬ 
ite  material ; 

t  t  t 

x'  y'  s  are  the  thicknesses  of  each  directional 
layer ; 

t  is  the  total  thickness  of  the  laminated 
plate; 

L,M,N  are  the  thickness  ratios  of  each  dir¬ 
ectional  layer; 

akj  is  the  stress  coefficient  (kj=x,y,s). 


Based  on  the  definitions  of  figure  1,  we  know  that  F  =F 

xx  a 

and  therefore  a  =1 . 

xx 


Fig.  2.  Schematic  Drawing  of  Ultimate  Strength  for  the 
Laminated  Plate  in  Each  Direction. 

In  the  same  way,  we  can  obtain: 

F, -  (a„L  +  a„\l  -r  a,,N )  F. 

F,  =  {d„L + a,,M  -  a.,N )  F. 


Based  on  the  definitions  of  figures  3(a)  and  (b)  we  know 


that  a  =1 ,  a  =a  ,  a  =a 

yy  yx  xy  xs  ys 


a  =a  ,  and  a  and  a 7  are  not 
sx  sy  ss  ss 


necessarily  equal  yet  to  simplify  we  consider  that  a  =a  7 

SS  SS. 


Fig.  3.  Schematic  Drawing  of  Ultimate  Strength  for  the 
Laminated  Plate  in  Each  Direction. 
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Written  in  a  matrix  form,  it  is: 


‘  F.- 

Aw 

IF.  * 

F, 

- 

a9* 

MF. 

.  F,  . 

.  A/l  ^#f  A/l  . 

.  NF.  . 

If  we  use  the  ultimate  load  of 


the  unit  width  then  it  is: 


‘  D.  * 

a>>  o..  a„  ' 

‘  I.F.  - 

D, 

= 

a»«  an  0r, 

UF. 

.D.  . 

.  a/»  On  . 

.  *>F.  . 

IV.  The  Method  of  Lamination  Optimum  Design  Method 
1 .  Method  and  Steps 

Because  it  is  necessary  that  the  determined  variables  are 

t  ,  t  ,  t  ,  the  weight  function  of  the  non-moment  plate  is 
x  y  s 

f(t  .  t  ,t  )=(t  +t  +t  )AY.  In  the  formula,  A  is  the  plate  surface 
x  y  s  x  y  s 

area  and  Y  is  material's  specific  gravity 


The  constraint  condition  is  a  certain  strength  criterion. 
For  example,  using_the  Norris  criteria, 

*  <'■•'»'•)*( Br)‘  +(t£-)'  -iuM#)’-'-0- 

are  all  functions  of  t 

S> 

are  related  to  t  , t 


In  the  formula. 


and  then  N  , N 


D  ,D  and  D 
x  y  s 

and  N. 


,t  and  t  , 
x  y  s 


and  t 


Applying  the  Lagrangian  multiplier  method: 


W 


f  /y,  i.) 

(/,+/,+/,MY+ A  ((*£-;)  +(-j£)  " 


Separately  letting  -0,  £  =°  and  ^s=°'  W®  Can  attain  the  fo1 
lowing  three  formulas : 


AyD* 

2\F. 


AyDI 
2 \F. 


ArD\ 

2A>. 


a„Nl  +  af.Ni ^  j  —  0.5iV»iVy^a,^  £)'  )+°y*  (  )  ) 

a„jV;+«„JV’(-g‘  J-0.5iV.iV,[a,^^-)  +  a„  (-§*-)*] 
■a„iV»  +  a^iVj(-g-)*-0.5^.Ary[a„(-g^)  +  ay,  (  g;)‘ ) 


+  a, 


(1) 

(2) 

(3) 


+a„.Vi( 


When  formula  (1)  is  reduced  into  formula  (2)  and  the  agx=agy 
relation  is  applied,  we  obtain: 

(a„-a„)jV;+  -0.5iV.N,[Ca..-a„)  (g;  )+<«■»-«»>  (f^)  ]"  0 


When  formula  (1)  is  reduced  into  formula  (3),  we  obtain: 


<°"  -O.SN.N,[(at,- +  (a,, _ «*)(-§“-)* J 


(5) 


The  selected  strength  criteria  can  be  rewritten  as: 


The  found  in  formula  (4)  is  substituted  into  formula  (5) 

7  /Dx\ 

and  then  the  found  is  substituted  into  formula  (6)  and  we 

s  of  Dx 

proceed  to  find  D  .  The  positive  and  negative  signs. are  identical 

X  A 

to  those  of  Nx. 

We  have  already  obtained  from  the  last  section  that: 


If, 


When  the  found  Dx,  —  and  are  substituted  into  formula 

y  s 

(7),  if  00  F  are  already  known,  then  we  can  find  the 

3 

values  of  t  , t  „  and  t  .  When  t.  is  divided  by  the  thickness  of 
x  y  s  k 

a  single  Layer  plate  then  we  obtain  the  layer  number  n^  of 
each  directional  layer.  However,  in  most  situations  the  n^  is 
not  nffcessai' ily  a  positive  integer  so  that  at  this  time  we  can 
take  the  closest  positive  integer  as  the  feasible  scheme  value. 
In  this  i?ay,  we  have  several  sets  of  feasible  schemes  and  its 
safety  tolerance  for  each  set  of  calculations  is 


The  selection  of  M.S.  is  the  largest  positive  value  as  the 
optimum  design  scheme. 


2.  Calculation  Examples 

A  cylindrical  structure  having  a  radius  of  R  *  10  inches,  a 
maximum  pressure  inside  the  cylinder  of  P»  1000  pounds/inch2  also 
sustains  a  maximum  torque  of  2  million  inch  pounds.  The  material 
used  is  boron/epoxy  composite  material. 


its  single  direction  permissable  strength  is  Fa=173,000  pounds/ 

2 

inch  ,  and  the  single  layer  thickness  of  each  type  of  lamination 
(0°,  90°,  i.45°)  is  tQ=0.0075  inches.  After  arranging  the  test 
data,  the  stress  coefficient  matrix  which  is 


On  0„ 

1  -0.1266  0.2133  ' 

a,,  a„  a,. 

- 

-0.1266^)*  1  0.2133 

.  0/»  Ott  . 

0.2133  0.2133  0.4266  . 

finds  the  lamination  optimum  design  scheme  of  this  structured 
cyl inder . 


Fig.  4.  Example  Figure 


The  solution  of  the  external  load  is 

Nx=PR=l  x  10=10,000  pounds/inch 

Ny=Pj=5 , 000  pounds/inch 
T 

N  = - --=3.183  pounds/inch 

s  27TR 


N 


N 

^=0.3183 


therefore 


1 

-0.1266 

0.2133 

'  1.202 

0.304 

-0.735  ' 

-0.0316 

1 

0.2133 

C  AY1  = 

0.181 

1.166 

-0.674 

.  0.2133 

0.2133 

0.4266 

. 

• 

.-0.674 

-0.735 

3.049  J 

When  already  known  data  is  substituted  into  formulas  (4) ,  (5) 
and  (6),  we  can  obtain  52=1.772,  5s=2.362  and  D  =12,100  pounds/ 
inch.  Further,  when  the  obtained  results  are  substituted  into 
formula  (7) ,  we  can  obtain 


//I  r  1.202  0.304  -0.735  If  1  1  TO  0743' 

!  |2,o  1 

■  1  1  0.181  1.166  -0.674  |  1.772  0.0387  Ifi  J‘ 

I  (1) 

-0.674  -0.735  3.049  A  2  362  J  LO.OUIJ 


1  I  (1) 

3.049  .L  2.362  J  LO.OUIJ 


9.9' 

5.2  J2  (2) 


f0.0743' 

if  =  0.0075  0  0387 
./,!  1.0.0141. 


Key:  1.  Inch 
2 .  Layer 


Because  2n=16.9  layer,  therefore  we  take  the  Zn=17  layer.  In 
this  way  we  have  the  following  three  types  of  feasible  schemes 
and  separately  calculate  their  safety  tolerances. 

D.  -  Ft  ( o,J,  +  <hj,- + a„l ,)  =*  1 . 2  9 8»i, -  0 . 1 6  +  0 . 27 7n, 

D,“  F.(a,.t, •+• a^, +a„/,)  -  -  0 . 04  In,  + 1 .  298«, + 0 . 277», 

D,  -  + a,J, + a„t,)  «  0 . 277  («, + n. + 2n.) 


(1)  «H-»* 


Based  on  the  M.S.  values,  we  should  select  scheme  2  as  the 
lamination  optimum  design  scheme. 

3 .  Improvement  of  Method 

During  actual  operation,  the  calculation  of  the  structural 
design  is  a  step-by-step  approximation  process.  In  addition,  the 
external  loads  and  stressed  panels  often  require  partitioned 
calculations  and  thus  in  order  to  decrease  the  amount  of  work 
and  make  it  convenient  for  using  an  electronic  computer,  we  can 
use  the  above  mentioned  calculation  measures  to  make  the 
following  changes.  According  to  the  matrix  formula  we  know 


A„ 

A., 


Aym  Atm 

An  Att 
Ayg  Ait 


In  the  formula,  a*a,ft  etc.  Formula 

(7)  can  be  rewritten  as: 


r#.l 

AnM  Afu  Atu 

r  a  1 

J 

Dm 

~FJa \ 

A»W  Afy  Atf 

D, 

D. 

n 

j.. 

Amg  Aft  Aft 

JL. 

-  D. 

Therefore,  we  obtain 


I, 

U 


Af+Ayy 

Ay  +  Aym-^  +  Aim-fy- 


^1/  ^  Aft  r\  ^  At 


An+Am 


A 

D. 

by 


D. 

Dm 

~D, 


Dm  +A,‘U. 


(8) 


(9) 
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Furthermore,  formulas  (4) ,  (5)  and  (6)  can  separately  be  re¬ 
written  as: 


(10) 

(fry<^  (&y) 

(zr)  - 0  (id 

“  - 1  +($/iI(frJ  <-fcUMk)W 


When  the  (ft  from  formula  (10)  is  substituted  into  formula 

(11)  and  the  ft).  of  formula  (11)  is  then  substituted  into 
formulas  (8)  ani  (9),  we  can  find  r1  and  r2-  At  the  same  time, 
when  ft)  and  | are  substituted  into  formula  (12)  ,  we  can 


when 


find 

Nx 


Letting  s*=R.  then  D=RN 


x  x 


Further 


D,  —  + atyi,+ a,,/,)  =•  Fj,(a„ + a„r, + a„rt) 


Therefore ,  we  obtain 


RN.--F  J.(a„ + a.»Ti + a  „rs) 


We  know  from  this  that 


FJ  ,  F.  (  .ntm_Ri±+n+rti 

Lm-ll-m - 1 - 

t  1  +r1+rI 

WsJi-s - Si - 

"  t  1  +r,+r. 
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Only  if  we  give  the  values  of  5/  and  **-,  based  on  the 

x  x 

known  (a)  and  the  F  in  accordance  with  formulas  (8) -(15),  we 

can  find  the  corresponding  numerical  values  of  L  and  N,  and 

Ft  N  N 

If  we  take  a  series  of  ^  and  ^  based  on  a  certian  specific 
N  N  N 

value  then  we  can  use  hand  computation  or  an  electronic  computer 
to  calculate  the  following  table  for  the  approximate  step-by- 
step  calculation  during  the  preliminary  design. 


4.  Special  Circumstances 

(1)  If  there  is  only  the  action  of 

The  optimum  design  of  lamination  is  naturally  accomplished 
by  the  lamination  along  the  k  direction. 


(2)  If  there  is  only  the  action  of  Nx  and  Ng  then  .V,°»  0 
This  type  of  loaded  situation  has  practical  significance  for  the 
stress  panel  of  the  wing  surface  and  therefore  we  will  discuss 
this  in  greater  detail. 


From  formula  (11)  we  can  obtain: 


(16) 


Because 


we  therefore  obtain 


(17) 


From  formula  (12)  we  can  obtain 


(18) 


From  formula  (13) 


Ft  R(i+r > 
N,  ~  a„+a..r 


(IS) 


From  formula  (14) 
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L 


1 

1  +  r 


(20) 


and  from  formula  (15) 


N 


r 

1  +  r 


1  ~L 


Formulas  (16) -(20)  are  very  easily  arranged  into  computer  pro¬ 
grams.  For  example,  if  the  stress  coefficient  and  F^  are  taken 
to  be  the  numerical  values  of  the  above  examples  after  going 
through  an  electronic  computer  we  can  obtain  the  following 
table 


■V,  j 

•  W7  ! 

0 

0.1 

o.z 

0.3 

0.4 

0.5 

0.6 

0.7 

1 

0.8  1 

| 

' 

0.9  1.0 

Fat 

«•  1 

<1) 

1.043 

t 

1.370 

1.667 

1.947 

2.217 

2.481 

!  2.738 

2.992 

3.242  \  3.489 

L  ;  (1) 

1.229 

0.979  0  926  0.718  0.636  0.371  '  0.51T 

1,11 

. 

0.472  1  0.433  i  0.399 

1  1  i 

If  the  data  in  the  table  are  drawn  into  figure  lines  then  it  is 
as  shown  in  fig.  5. 
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Fig.  5.  The  N  /N  -F  t/N  -L  Relation  Curve 


N  D 

When  g^  is  very  small,  we  know  from  formula  (16)  that  ^ 

X  N  D  s 

will  be  very  large  and  especially  when  g^=0,  will  become 

N  X  s 

infinite.  In  reality,  however,  when  rr^=0,  that  is  it  becomes  the 

Nx 

x  Ft 

single  load  of  N  =N  =0  and  N  ^0  then  there  should  be  and 

y  s  y  x  N 

N 

L=l.  Therefore,  the  L-  g^  curve  in  fig.  5  should  take  a  value  in 

X  N 

accordance  with  the  dotted  line  when  ^  is  smaller  than  a  certain 


value  (for  example,  taking  g^=0.3).  The  dotted  line  area  shows 


the  transformations  within  two  types  (single  load  and  double  load) 
of  minimum  weight  designs.  Because  of  the  selection  of  the  dotted 
line  area's  curve  form  and  stopping  point  there  is  a  certain 
randomness  and  therefore  its  error  can  be  relatively  large. 
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V.  Discussion  of  Technique 

1.  Can  we  differentiate  and  determine  the  number  of  layers 
of  each  directional  layer  based  on  the  external  load  of  each 
direction? 

Also,  is  it  possible  to  carry  out  a  lamination  design  based 
Nk 

on  the  t  =—— w-  (k=x ,  y ,  s)  formula.  When  we  still  take  the  above 
k  akk  a 

example's  data,  we  can  obtain: 


<,s 


^-•.0578 AM. 


(2) 


0.0075 


U 

i. 


5 

173 


0.0289 


0.0289 

0.0075 


3.183 

0.4266X173 


0.043^-t, 


n. 


3.85-  4 

0.0431 _ 
0.0075“ 


m (4) 

(6) 

5.74—  6  JS 


Key:  1.  Inch 

2 .  Layer 

3 .  Inch 

4 .  Layer 

5 .  Inch 

6 .  Layer 


We  know  that  the  differences  of  the  optimum  schemes  in  the 
contrasting  examples  are  relatively  large  and  it  appears  that 
the  main  reason  for  these  differences  is,  as  regards  the 
laminated  plate,  causing  the  external  load  to  be  single  direction 
al.  The  stress  conditions  of  each  directional  layer  are  still 
complex  stress  conditions  (if  we  do  not  calculate  the  stress 
in  the  layer  then  they  are  plane  stress  conditions)  and  the 
interaction  within  each  directional  layer  cannot  be  overlooked. 


2.  Stress  coefficient  a.,  and  single  direction  composite 


longitudinal  tensile  strength  F 

We  can  see  from  the  numerical  examples  that  the  influence 
of  the  numerical  value  of  a.  .  on  the  calculation  results  is 

quite  large.  The  value  of  a^j  is  obtained  from  special  single 

load  tests  based  on  the  different  percentages  of  each  layer  of 
the  composite  on  the  x,  y  and  s  axes.  In  the  same  way,  F  is 
also  determined  from  the  tests. 

3 .  Strength  Data 

We  must  take  into  account  that  this  paper  used  the  laminated 

plate's  strength  criteria  and  that  this  was  for  a  non-singular 

laminated  plate.  Therefore,  each  of  its  data  items  should  use 

the  laminated  plate.  The  above  mentioned  formulas  were  based  on 

the  Norris  criteria  yet  if  we  use  the  Hill-Tsai  criteria  then 

the  formulas  can  have  fewer  differences  and  yet  its  results  can 

also  have  some  dissimilarities.  However,  the  two  formulas  are 

not  dissimilar  for  the  and  ^^0  situations.  If  we  use  the 

Nx  Nz 

Tsai-Wu  tensor  multinomial  criteria,  the  derived  formula  is 
quite  complex  and  inconvenient  to  use  because  there  are  very 
many  number  items  and  also  because  it  is  necessary  to  use  the 
tensile  and  compression  limit  strength.  Therefore,  it  is  more 
convenient  to  use  the  previous  two  types  of  criteria  during  the 
preliminary  design. 
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Abstract 

When  the  upper  and/or  lower  panels  of  the  loading  box  of  the  aircraft 
wing  or  tail  are  made  of  the  fiber  reinforced  composite  laminates,  they  can 
frequently  be  simplified  as  a  non-moment  plate.  This  paper  introduces  an 
optimum  design  (i.  e.  minimum  weight  design)  procedure  of  the  laminated  plate 
on  static  failure  strength  condition.  The  mathematical  tool  used  in  the  proce¬ 
dure  is  Lagrangian  Multiplier  method,  and  the  static  failure  strength  condition 
is  adopted  as  Hill~Tsai  criteria  or  Norris  criteria. 

The  formulae  for  optimum  design  have  been  not  only  derived,  but  also 
reformed  to  be  convenient  for  the  computer.  How  to  establish  the  ultimate 
strength  of  laminates  is  discussed  in  detail.  An  example,  illustrating  the 
solution  procedure  and  how  to  select  the  optimum  scheme  of  lamination  design, 
is  presented  in  the  paper.  Some  technical  problems  are  briefly  discussed  in  the 
last  part. 

At  the  stage  of  the  preliminary  structural  design,  this  procedure  can  be 
considered  as  an  engineering  method  of  lamination  optimum  design  for  the  loa¬ 
ding  panel  of  laminated  composite,  which  works  under  tension  or  compression 
(assuming  that  the  buckling  failure  would  not  occur). 
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